Boundary Detection using CWT - Interactive.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "03QpgZU3W2dg"
   },
   "source": [
    "# **Lithology Boundary Definition using CWT**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "nnev6gD1gXbX"
   },
   "source": [
    "Github : https://github.com/galena100/Transform2020\n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/galena100/Transform2020/master/t20-litho_boundary_from_gamma/test_on_synthetic_data.JPG\" alt=\"Your image title\" width=\"500\"/>\n",
    "\n",
    "We created this image from a synthetic GR log as part of this project.  It show the potential for finding lithology boundaries from geophysical logs.\n",
    "\n",
    " \n",
    "\n",
    "Based on \"***Improving Automated Geological Logging of Drill Holes by Incorporating Multiscale Spatial Methods***\" works by [Evelyn Jun Hill](https://www.researchgate.net/publication/339871641_Improving_Automated_Geological_Logging_of_Drill_Holes_by_Incorporating_Multiscale_Spatial_Methods). \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "xqyvrylOXAxy"
   },
   "source": [
    "**1. What did you guys do?**\n",
    "---\n",
    "Estimate lithology boundaries from a downhole geophysical log, using a continuous wavelet transform (CWT).\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "mJB7MwtfYZZ5"
   },
   "source": [
    "## **2. Who are you guys?**\n",
    "\n",
    "- **Jared Armstrong** :\n",
    " \n",
    "  \"*I love working at the intersection of mining coal geology and Database Technology.  Free your data to dissolve your difficulties*\"\n",
    "\n",
    "  E-mail : jared@wholesolutions.com.au\n",
    "\n",
    "  Slack : @Jared \n",
    "\n",
    "- **Leo C. Dinendra** : \n",
    "\n",
    "  \"*A geophysicist DM worker in an oil company. I love for-loops and will fight for it.*\"\n",
    "\n",
    "  E-mail : leocd91@gmail.com , [LinkedIn](https://www.linkedin.com/in/leo-c-0988727b/)\n",
    "\n",
    "  Slack : @leo \n",
    "\n",
    "- **Martin Bentley** :\n",
    "\n",
    "  \"*Messing about with computers and rocks. Accidental geophysicist.*\"\n",
    "\n",
    "  E-mail : astonmartin.bentley@gmail.com\n",
    "\n",
    "  Slack : @mtb-za\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "UwX92F4YYGBX"
   },
   "source": [
    "**3. How you guys do it?**\n",
    "---\n",
    "\n",
    "Run all of the cells until the **Making Things Happen** heading.\n",
    "\n",
    "### Standard Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "sG1uV6iAWI_7"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.signal import savgol_filter\n",
    "import pywt\n",
    "import ipywidgets as widgets\n",
    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
    "from IPython.display import display"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining our functions\n",
    "\n",
    "### Importing data and getting our lithology set up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "oOyB2mIUDWK8"
   },
   "outputs": [],
   "source": [
    "def extract_data(fname, borehole_id='test001', data_type='DENB_g_cc'):\n",
    "    '''\n",
    "    Read a CSV file. Extract only data for `borehole` in dataset, then use only `data_type` for further plots.\n",
    "    returns:\n",
    "        data: series of interest\n",
    "        curve: Sav-Gol filtered dataset\n",
    "    '''\n",
    "    raw_data = pd.read_csv(fname)\n",
    "    bore = raw_data.loc[raw_data['HOLEID'] == borehole_id].copy() # get data for the requested hole.\n",
    "    bore.dropna(subset=[data_type], inplace=True)\n",
    "    bore['DEPTH'] = bore['DEPTH'] * -1 #  Make depths positive down.\n",
    "    bore['savgol'] = savgol_filter(bore[data_type], 21, 2)                                      #getting the GR data and filter using savgol\n",
    "    \n",
    "    return bore, borehole_id, data_type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We need a dict to convert lithologies to colours.\n",
    "lith_colours = {\n",
    "    'Basalt': 'purple',\n",
    "    'Calcite': 'lightgrey',\n",
    "    'Carbonaceous Mudstone': 'lightgreen',\n",
    "    'Clay': 'green',\n",
    "    'Claystone': 'darkgreen',\n",
    "    'Coal': 'black',\n",
    "    'Core Loss': 'white',\n",
    "    'Mudstone': 'darkgrey',\n",
    "    'Sand': 'gold',\n",
    "    'Sandstone': 'orange',\n",
    "    'Siderite': 'lightblue',\n",
    "    'Siltstone': 'brown',\n",
    "    'Soil': 'chocolate',\n",
    "    'Tuff': 'yellow',\n",
    "    'Intermediate Intrusive': 'rebeccapurple',\n",
    "    'Pyrite': 'darkgoldenrod',\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_lithology(borehole='test001', lithcolours=lith_colours):\n",
    "    '''\n",
    "    Using the provided lithology and lookup table, will give back a lithology table\n",
    "    and list of colours.\n",
    "    \n",
    "    Needs a dictionary of lithologies with associated colours.\n",
    "    '''\n",
    "    lithology = pd.read_csv('t20-litho_boundary_from_gamma/2_litho.csv')\n",
    "    lithology_lookup = pd.read_csv('t20-litho_boundary_from_gamma/2_Lookup_Table.csv')\n",
    "    descriptions = lithology_lookup.to_dict()['DESCRIPTION'].values()\n",
    "    lookups = lithology_lookup.to_dict()['LOOKUP'].values()\n",
    "\n",
    "    lithology_lookup = {}\n",
    "    for lookup, description in zip(lookups, descriptions):\n",
    "        lithology_lookup[lookup] = description\n",
    "    lithology.replace(to_replace=lithology_lookup, value=None, inplace=True)\n",
    "    litho_first_hole = lithology[lithology['HOLEID'] == borehole].copy()\n",
    "    litho_first_hole['spans'] = litho_first_hole['GEOLTO'] - litho_first_hole['GEOLFROM']\n",
    "    \n",
    "    litho_first_hole['colour'] = litho_first_hole['VALUE']\n",
    "    litho_first_hole['colour'].replace(lith_colours, value=None, inplace=True)\n",
    "    \n",
    "    colours = litho_first_hole['colour'].to_list()\n",
    "    \n",
    "    return litho_first_hole, colours"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "4i2wHF8l7A0h"
   },
   "source": [
    "### CWT, Boundary Definition, and Boundary Strength Computation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### These are support functions, we need them for the data we plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def resample_curve(curve, step):\n",
    "    '''Returns a pandas series, keeping every `step`th item.'''\n",
    "    return curve.iloc[list(range(0, len(curve), step))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_wavelet(curve, scale, wavelet_type, dt):\n",
    "    '''\n",
    "    Calculate the Continuous Wavelet Transform for a given CWT\n",
    "    Returns:\n",
    "        Transposed coefficients of the CWT\n",
    "    '''\n",
    "    [cfs, frequencies] = pywt.cwt(curve, scale, wavelet_type, dt)\n",
    "    return cfs.T, frequencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_cwt(curve, dt, fscales=np.arange(1, 100,1)):\n",
    "    '''Return coefficients of first and second order Gaussian CWT of given curve'''\n",
    "    zcwt2, g2_freq = make_wavelet(curve, fscales, 'gaus2', dt)\n",
    "    zcwt1, g1_freq = make_wavelet(curve, fscales, 'gaus1', dt)\n",
    "    \n",
    "    return zcwt2, zcwt1, g2_freq, g1_freq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_boundary_strength(cs, boundary_threshold=30):\n",
    "    '''\n",
    "    Generates a horizontal histogram of depths where the CWT is equal to zero. \n",
    "    '''\n",
    "    cs_allsegs = np.array(cs.allsegs[0])\n",
    "\n",
    "    #print(cs_allsegs[0][:])\n",
    "\n",
    "    #'''\n",
    "    y_counts = []\n",
    "    for cs_segs in cs_allsegs:\n",
    "        if cs_segs.shape[0] > boundary_threshold: #  We can use this for additional filtering to find boundaries at different wavelengths.\n",
    "            for seg in cs_segs:\n",
    "                y_counts.append(seg[1])\n",
    "    return y_counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 513
    },
    "colab_type": "code",
    "id": "9SokReraqXvk",
    "outputId": "a2dc5348-89c0-4ae1-8910-781b69c61fbb"
   },
   "outputs": [],
   "source": [
    "def plotting_function(borehole, plot_type, lithology,\n",
    "                      boundary_threshold=30,\n",
    "                      fscales=np.arange(1, 100,1),\n",
    "                      resample_curves=True,\n",
    "                      y_top=0, y_bot=100,\n",
    "                      borehole_id='A Borehole'):\n",
    "    # Useful variables:\n",
    "    dep = borehole[borehole['DEPTH'] > y_top]\n",
    "    dep = dep[borehole['DEPTH'] < y_bot]['DEPTH']\n",
    "    curve = borehole[borehole['DEPTH'] > y_top]\n",
    "    curve = curve[borehole['DEPTH'] < y_bot]['savgol']\n",
    "    plot_range = [y_top, y_bot]\n",
    "    \n",
    "    if resample_curves:\n",
    "        dep_resampled = resample_curve(dep, 10)\n",
    "        curve_resampled = resample_curve(curve, 10)\n",
    "    else:\n",
    "        dep_resampled = dep\n",
    "        curve_resampled = curve\n",
    "\n",
    "    dt = dep_resampled.iloc[1] - dep_resampled.iloc[0]\n",
    "\n",
    "    zcwt, zcwt2, freq1, freq2 = calculate_cwt(curve_resampled, dt)\n",
    "    \n",
    "    print('Making plots.')\n",
    "    # Begin plot\n",
    "    fig, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4, figsize=(12, 10))\n",
    "\n",
    "    # Data for each plot\n",
    "    # Base curve\n",
    "    ax1.plot(curve_resampled, dep_resampled)\n",
    "    # CWT\n",
    "    cs=ax2.contour(fscales, dep_resampled, zcwt, levels=[0.0], colors='k', extend='both')\n",
    "    print('Calculating boundary strengths.')\n",
    "    \n",
    "    ## We need to calculate our boundary strengths\n",
    "    y_counts = make_boundary_strength(cs, boundary_threshold)\n",
    "    \n",
    "    ## Now carry on plotting.\n",
    "    ax2.contourf(fscales, dep_resampled, zcwt,\n",
    "                 levels=np.arange(np.min(zcwt), np.max(zcwt), 0.5), \n",
    "                 cmap=\"RdBu_r\", extend='both')\n",
    "    # Boundary strength\n",
    "    ax3.hist(np.array(y_counts), bins=int(np.floor(dep.max())), orientation='horizontal', color='k')\n",
    "    # Lithology\n",
    "    ax4.barh(lithology['GEOLFROM'], align='edge',\n",
    "             height=list(lithology['spans']),\n",
    "             width=10, color=lithology['colour'],)\n",
    "\n",
    "    # Set plot limits\n",
    "    ax1.set_ylim(plot_range)\n",
    "    ax2.set_ylim(plot_range)\n",
    "    ax3.set_ylim(plot_range)\n",
    "    ax4.set_ylim(plot_range)\n",
    "    ax4.set_xlim([0.0,1.0])\n",
    "\n",
    "    # Invert axes\n",
    "    ax1.invert_yaxis()\n",
    "    ax2.invert_yaxis()\n",
    "    ax3.invert_yaxis()\n",
    "    ax4.invert_yaxis()\n",
    "\n",
    "    # Title\n",
    "    ax1.set_xlabel(f'Smoothed {plot_type}')\n",
    "    ax2.set_title(f'CWT from {np.round(np.min(freq1),2)} Hz to {np.round(np.max(freq2),2)} Hz')\n",
    "    ax3.set_title('Boundary Strength')\n",
    "    ax4.set_title('Lithology')\n",
    "\n",
    "    # Axis Labels\n",
    "    ax1.set_ylabel('Depth (m)')\n",
    "    ax2.xaxis.set_visible(False)\n",
    "    ax2.yaxis.set_visible(False)\n",
    "    ax3.xaxis.set_visible(False)\n",
    "    ax3.yaxis.set_visible(False)\n",
    "    ax3.set_xlabel('Boundary Strength')\n",
    "    ax4.xaxis.set_visible(False)\n",
    "    #ax4.yaxis.set_visible(False) #  Keeping this visible, because the spine length goes weird sometimes.\n",
    "\n",
    "    ax4.set_aspect(0.05)\n",
    "    \n",
    "    plt.title(f\"CWT of {plot_type} in {borehole_id}\")\n",
    "    \n",
    "    # Adding gridlines\n",
    "    ax1.grid(zorder=1)\n",
    "    ax2.grid(zorder=4);ax2.set_axisbelow(True)\n",
    "    ax3.grid(zorder=4);ax3.set_axisbelow(True)\n",
    "    ax4.grid(zorder=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Making things happen:\n",
    "\n",
    "Run the next cell, and choose your data. There are three boreholes, and four well curves to choose from."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fname = 't20-litho_boundary_from_gamma/2_geophys.csv'\n",
    "method = 'MC2F_us_f'\n",
    "borehole_id = 'test001'\n",
    "select_data = interactive(extract_data,\n",
    "                          borehole_id=['test001', 'test002', 'test003'],\n",
    "                          data_type=['GRDE_gapi', 'DENB_g_cc', 'MC2F_us_f', 'CADE_mm'],\n",
    "                          fname=fixed(fname))\n",
    "display(select_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that you have chosen your data, set up some of the plot's parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "borehole = select_data.result[0]\n",
    "\n",
    "y_top_w = widgets.IntSlider(value=int(borehole['DEPTH'].min()),\n",
    "                            min=int(borehole['DEPTH'].min()),\n",
    "                            max=int(borehole['DEPTH'].max()-5),\n",
    "                            step=5, description='Window Top')\n",
    "y_bot_w = widgets.IntSlider(value=int(borehole['DEPTH'].max()-5),\n",
    "                            min=5, max=int(borehole['DEPTH'].max()),\n",
    "                            step=5, description='Window Bottom')\n",
    "threshold_slider = widgets.IntSlider(value=30, min=0, max=100, step=5, description='Boundary Threshold')\n",
    "\n",
    "y_top_w_link = widgets.dlink((y_top_w, 'value'), (y_bot_w, 'min'))\n",
    "y_bot_w_link = widgets.dlink((y_bot_w, 'value'), (y_top_w, 'max'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can tweak the sensitivity of the edge detection here, by changing the 'Boundary Threshold'. The 'Window Top' and 'Window Bottom' will let you zoom in on sections of the borehole. Sometimes the lithology plot goes a little screwy when doing this. No idea why.\n",
    "\n",
    "If you want a different curve, or a different borehole, go back to the selection cell, then run the cell above and below this text again. Have fun!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "lithology, colours = make_lithology(borehole=borehole_id)\n",
    "plot_data = interact_manual(plotting_function,\n",
    "                            borehole=fixed(select_data.result[0]), \n",
    "                            plot_type=fixed(select_data.result[2]), \n",
    "                            lithology=fixed(lithology),\n",
    "                            boundary_threshold=threshold_slider,\n",
    "                            fscales=fixed(np.arange(1, 100,1)),\n",
    "                            resample_curves=True,\n",
    "                            y_top=y_top_w,\n",
    "                            y_bot=y_bot_w,\n",
    "                            borehole_id=fixed(borehole['HOLEID'].iloc[0]),\n",
    "                            )"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "T20-Litho_boundary_from_GR_use_CWT.ipynb",
   "provenance": []
  },
  "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}