Getting_started_with_Python.ipynb

{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["# Tutorial: Getting Started with Python\n", "\n", "DATE: Monday June 8, 2020 08:00 - 11:00 UTC\n", "\n", "AUDIENCE: Beginner\n", "\n", "INSTRUCTOR: Robert Leckenby, Digital Geoscientist, [Agile](https://agilescientific.com/)\n", "\n", "### Video Stream:\n", "<br>\n", "<a href=\"https://youtu.be/iIOMiN8Cacs\" target=\"_blank\">\n", "    <img src=\"./youtube_thumbnail.png\" width=800 />\n", "</a>\n", "\n", "## Welcome\n", "Welcome to Transform 2020 and to Getting Started with Python! This session is open to all but is aimed at people who have never coded at all or who don't know Python. Ideally you should have watched the setup videos for [Windows](https://youtu.be/FdatS_NKVrM) or [Linux](https://youtu.be/3ncwbHyZeAg) and/or followed the instructions in the [Installation guide for the tutorials](http://swu.ng/t20-python-setup) that contained the information in those videos and more. This session is meant as a demo but if you followed those installation steps, you should be able to follow along. \n", "\n", "All session details are available [here](https://transform2020.sched.com/).\n", "\n", "## Agenda\n", "START:&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;08.00 UTC\n", "\n", "BREAK:&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;09.30 UTC\n", "\n", "RESUME:&nbsp;&nbsp;&nbsp;09.40 UTC\n", "\n", "END:&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;11.00 UTC\n", "\n", "## Content\n", "\n", "- Simple math, logic and flow control\n", "- Modelling [Darcy's law](https://en.wikipedia.org/wiki/Darcy%27s_law) with a function\n", "- Computing fold wavelengths with the [Ramberg-Biot](https://blogs.agu.org/mountainbeltway/2010/10/15/friday-fold-wavelength-contrast/) equation\n", "- Geophysical wavelets with [`bruges`](https://github.com/agile-geoscience/bruges)\n", "- Mapping DEM's with [`matplotlib`](https://matplotlib.org/)\n", "- Data wrangling and QC well data with [`welly`](https://github.com/agile-geoscience/welly)\n", "- Scanning a seismic cube with [`segyio`](https://github.com/equinor/segyio)\n", "\n", "## Set up\n", "\n", "The [anaconda documentation](https://docs.conda.io/projects/conda/en/latest/user-guide/concepts/environments.html#conda-environments) tells us that: \"A conda environment [...] contains a specific collection of conda packages that you have installed.\" We use these environments to isolate different projects.\n", "\n", "We will be using an environment called `transform2020` (but you can call it whatever you wish) that was created with:\n", "\n", "`conda create --name transform2020 python=3.7 anaconda`\n", "\n", "We then also install the [welly](https://github.com/agile-geoscience/welly), [segyio](https://github.com/equinor/segyio) and [bruges](https://github.com/agile-geoscience/bruges) packages to import \\*.las and \\*.sgy files and to perform geophysical calculations. First we activated the environment:\n", "\n", "`conda activate transform2020`\n", "\n", "and then installed the additional libraries:\n", "\n", "`pip install welly segyio bruges`\n", "\n", "You can then access this environment inside a jupyter notebook using one of two methods:\n", "1. By running `pip install nb_conda_kernels` in `base` and then running `jupyter notebook` with `base` activated.\n", "2. By running `python -m ipykernel install --user --name transform2020` after creating the `transform2020` environment and then running `jupyter notebook` with either `base` or `transform2020` activated.\n", "\n", "In both cases you then need to select the `transform2020` environment by selecting `Kernel>Change kernel>kernel name` in the jupyter notebook browser window."]}, {"cell_type": "markdown", "metadata": {}, "source": ["## Python basics: numbers and functions"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "markdown", "metadata": {}, "source": ["Let's express an equation symbolically, then in Python.\n", "\n", "$$ t = xy + z $$"]}, {"cell_type": "markdown", "metadata": {}, "source": ["<div style=\"padding:0.5em 1em 1em; background-color:#eeffee; border-radius:5px; border:1px solid #77aa77;\">\n", "\n", "<h2>EXERCISE</h2>\n", "\n", "<p><b>Write the equation as a line of Python code, using the same symbols.</b></p>\n", "\n", "<p><b>You'll also need to define some variables. Use 3, 10, and 1 for `x`, `y`, and `z` respectively.</b></p>\n", "\n", "</div>"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Here's how that looks in Python:"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "markdown", "metadata": {}, "source": ["Fine, but it's not very convenient to have to define `x`, `y`, and `z`, then write out the equation. If I want to do it again on different values of `x`, `y`, and `z`, then I'll have to write it all out again."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["x = 4\n", "y = 100\n", "z = 0.1\n", "\n", "t = x * y + z\n", "\n", "t"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Instead, we can define a function:"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["compute_t(3, 10, 1)"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["compute_t(4, 100, 0.1)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["There's no reason we can't pass variables into here:"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["compute_t(a, b, 5)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["What if we want to compute a range of values for `x`?"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["compute_t(x, 10, 1)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Pretty convenient!\n", "\n", "And what if we need to check for conditions?\n", "\n", "We can use `if ... else` statements for this:"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "markdown", "metadata": {}, "source": ["## Transforming geoscience equations into Python code\n", "\n", "We will take some classical geoscience equations and see how to code them in Python, both for single-use and in `functions` for reuse."]}, {"cell_type": "markdown", "metadata": {}, "source": ["### Darcy's law\n", "\n", "[Darcy's law](https://en.wikipedia.org/wiki/Darcy%27s_law) describes flow rate in porous media:\n", "\n", "$$ Q = \\frac{kA}{\\mu L} \\Delta P$$\n", "\n", "where:\n", "\n", "- $Q$ = flow rate [m<sup>3</sup>/s]\n", "- $k$ = permeability [m<sup>2</sup>]\n", "- $A$ = cross-sectional area [m<sup>2</sup>]\n", "- $\\mu$ = viscosity [Pa.s]\n", "- L = length of sample [m]\n", "- $\\Delta P$ = pressure difference accross the sample [Pa]\n", "\n", "All units are SI."]}, {"cell_type": "markdown", "metadata": {}, "source": ["<div style=\"padding:0.5em 1em 1em; background-color:#eeffee; border-radius:5px; border:1px solid #77aa77;\">\n", "\n", "<h2>EXERCISE</h2>\n", "\n", "<p><b>Implement Darcy's law in Python. You can do it either as a 'raw' Python expression, like the `t = x * y + z` example above, or as a function, like the `compute_t()` example.</b></p>\n", "</div>"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "markdown", "metadata": {}, "source": ["Alternatively, using a Python function:"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "markdown", "metadata": {}, "source": ["Let's use some more realistic values.\n", "\n", "Let's consider a side wall core:\n", "\n", "End view | Side view\n", ":-- | :-- \n", "![alt](./data/swc_end.png) | ![alt](./data/swc_side.png)"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["k = 4.9346165e-13 # ~500 mD\n", "A = 2.03e-3       # 1 inch diameter plug\n", "mu = 8.9e-4       # ~dynamic viscosity of water at about 25\u00b0C\n", "L = 0.0508        # 2 inch plug\n", "dP = 689476       # ~100 psi\n", "\n", "Q = flow_rate(k, A, mu, L, dP)\n", "\n", "print(f'The flow rate `Q` is approx: {Q:.1e} m3/s or {Q * 3600:.2f} m3/hr')"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Let's use a NumPy array to make a plot of flow rate vs plug length. We'll need to define the range of plug lengths:"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["Q = flow_rate(k, A, mu, L, dP) * 3600"]}, {"cell_type": "code", "execution_count": null, "metadata": {"scrolled": false}, "outputs": [], "source": ["%matplotlib inline\n", "import matplotlib.pyplot as plt  \n", "\n", "plt.plot(L, Q)\n", "plt.show()"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["from ipywidgets import interact\n", "\n", "@interact(dP=(0, 4, 0.1))\n", "def make_plot(dP):\n", "    \n", "    # We have to compute Q every time the widget changes...\n", "    Q = flow_rate(k, A, mu, L, dP*1e6) * 3600\n", "    \n", "    # ...and then make the plot again.\n", "    fig, ax = plt.subplots(figsize=(16,6))\n", "    ax.plot(L, Q)\n", "    ax.grid()\n", "    ax.set_title('Flow rate [m3/hr] versus plug length [m]', fontsize=16)\n", "    ax.set_ylabel('Flow rate [m3/hr]', fontsize=14)\n", "    ax.set_xlabel('Plug length [m]', fontsize=14)\n", "    ax.text(0.83, 0.088, f\"dP: {dP} MPa\", size=20)\n", "    ax.set_ylim(0, 0.1)\n", "    plt.show()\n", "    \n", "    return"]}, {"cell_type": "markdown", "metadata": {}, "source": ["### Ramberg-Biot Equation\n", "\n", "The [Ramberg-Biot equation](http://www.files.ethz.ch/structuralgeology/JPB/files/English/7folding.pdf) describes the wavelength $L$ of a fold:\n", "\n", "$$ L = 2 \\pi T \\left(\\frac{\\nu}{6\\nu_0}\\right)^{1/3}$$\n", "\n", "where:\n", "\n", "- $L$ = wavelength of the fold [m]\n", "- $\\pi$ = 3.14 [unitless]\n", "- $T$ = folded layer thickness [m]\n", "- $\\nu$ = viscosity of folded layer [Pa.s]\n", "- $\\nu_0$ = viscosity of cross-cut layer [Pa.s]\n", "\n", "An [example](https://blogs.agu.org/mountainbeltway/2010/10/15/friday-fold-wavelength-contrast/) of the Ramberg-Biot equation in action.\n", "\n", "This time we will define the function directly, notice we are using the $\\pi$ constant from `numpy`.\n", "\n", "First an example to discuss what we are looking into:"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["fig, ax = plt.subplots(figsize=(16,4))\n", "ax.plot(-1*np.sin(np.sin(np.linspace(0,3*np.pi,100))),lw=4, c='g')\n", "ax.plot(-1*np.sin(np.linspace(0,5*np.pi,100)) - 2,lw=4, c='b')\n", "ax.set_yticklabels([])\n", "ax.set_xticks(np.arange(0,110,10))\n", "ax.grid()\n", "\n", "ax.set_title('Simplified illustration of Ramberg-Biot usage', fontsize=16)\n", "textstr_0 = 'Long wavelength in high\\nviscosity enclosing medium'\n", "textstr_1 = 'Short wavelength in low\\nviscosity enclosing medium'\n", "props_cret = dict(boxstyle='square', facecolor='lightgreen', alpha=0.4)\n", "props_jur = dict(boxstyle='square', facecolor='lightblue', alpha=0.4)\n", "props_well = dict(boxstyle='square', facecolor='lightgrey', alpha=0.4)\n", "\n", "ax.axvline(50, lw=3, c='k', ls=':')\n", "\n", "ax.text(0.09, 0.95, textstr_0, transform=ax.transAxes, fontsize=14,\n", "        verticalalignment='top', bbox=props_cret)\n", "ax.text(0.01, 0.48, textstr_1, transform=ax.transAxes, fontsize=14,\n", "        verticalalignment='top', bbox=props_jur)\n", "ax.text(0.7, 0.78, 'Cretaceous', transform=ax.transAxes, fontsize=14,\n", "        verticalalignment='top', bbox=props_cret)\n", "ax.text(0.77, 0.4, 'Jurassic', transform=ax.transAxes, fontsize=14,\n", "        verticalalignment='top', bbox=props_jur)\n", "ax.text(0.48, 0.55, 'Well', transform=ax.transAxes, fontsize=14,\n", "        verticalalignment='top', bbox=props_well, rotation='vertical')\n", "\n", "plt.show()"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Gara Anticline - Kurdistan (satellite) | Gara Anticline - Kurdistan (terrain)\n", ":-- | :-- \n", "![alt](./data/gara_satellite.png) | ![alt](./data/gara_terrain.png)\n", "\n", "[Google map source location](https://www.google.ch/maps/@36.9742638,43.6135018,12z/data=!5m1!1e4)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["<div style=\"padding:0.5em 1em 1em; background-color:#eeffee; border-radius:5px; border:1px solid #77aa77;\">\n", "\n", "<h2>EXERCISE</h2>\n", "\n", "<p><b>Write the Ramberg-Biot equation in a Python function, using the same variable names as shown above.</b></p>\n", "\n", "\n", "\n", "</div>"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "markdown", "metadata": {}, "source": ["Given an input data set with the required arguments:\n", "\n", "|Period    |Formation                   |Thickness [m]|$\\nu$ [Pa.s]|$\\nu_0$ [Pa.s]|\n", "|:---------|:---------------------------|:-------:|:---:|:-----:|\n", "|Cretaceous|Aqra Bekhme Kometan Qamchuqa|730      |19   |17     |\n", "|Jurassic  |Naokelekan Sargelu          |80       |17   |19     |\n", "|Jurassic  |Sekhanian Sargelu           |420      |19   |17     |\n", "|Triassic  |Baluti                      |30       |17   |19     |\n", "|Triassic  |Kurra Chine                 |834      |18   |17     |\n", "|Triassic  |Geli                        |575      |18   |17     |\n", "|Triassic  |Beduh                       |64       |17   |19     |\n", "|Triassic  |Chia Zara                   |811      |19   |17     |\n", "\n", "***Nota Bene***: The thickness values were measured from wells, the viscosity values are simple examples values but are not measured.\n", "\n", "We could print out the results of the calculation for each row, but we'll first need to store this table in a python data structure.\n", "\n", "We'll put the values together in a [`tuple`](https://docs.python.org/3.7/library/stdtypes.html?highlight=tuple#tuple), this will allow us to `index` into them to retrieve one value at a time. But then we'll put all of these values together in a dictionnary (a python [`dict`](https://docs.python.org/3.7/tutorial/datastructures.html#dictionaries)) where the keys will be the formation names:"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["data = {'ABKomQam':(730,19,17),\n", "        'NaokSargelu': (80,17,19),\n", "        'SekhSargelu': (420,19,17),\n", "        'Baluti': (30,17,19),\n", "        'KurraChine': (834,18,17),\n", "        'Geli': (575,18,17),\n", "        'Beduh': (64,17,19),\n", "        'ChiaZara': (811,19,17),\n", "       }"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "markdown", "metadata": {}, "source": ["Or save the results directly to a CSV file using the Python inbuilt [CSV module](https://docs.python.org/3.7/library/csv.html#module-csv):"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["import csv\n", "\n", "wavelengths = {}\n", "for formation, inputs in data.items():\n", "    T, nu, nu_zero = inputs\n", "    rb = ramberg_biot(T, nu, nu_zero)\n", "    wavelengths[formation] = rb\n", "\n", "with open('./data/wavelengths.csv', 'w') as f:\n", "    w = csv.writer(f)\n", "    w.writerows(wavelengths.items())"]}, {"cell_type": "markdown", "metadata": {}, "source": ["We could also have used [pandas](https://pandas.pydata.org/) to write _to_ disk but let's now read _from_ disk using it's `.read_csv` method:"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["import pandas as pd\n", "\n", "df_wave = pd.read_csv('./data/wavelengths.csv', names=['Formation', 'Wavelength'], index_col='Formation')\n", "df_wave"]}, {"cell_type": "markdown", "metadata": {}, "source": ["#### We could also plot the results:\n", "\n", "We can make a simple - but not very legible - plot with one line of code:"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["plt.bar(df_wave.index, df_wave['Wavelength'])"]}, {"cell_type": "markdown", "metadata": {}, "source": ["But now let's dress that plot up to make it more useful, first we'll set up a [`list`](https://docs.python.org/3.7/library/stdtypes.html?highlight=list#lists) to store [hex color codes](https://www.color-hex.com/), and then change the figure size, add a title, labels and a grid."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["colors = ['#339966', # Cretaceous\n", "          '#333399', # Jurassic\n", "          '#333399', # Jurassic\n", "          '#CC99FF', # Upper Triassic\n", "          '#CC99FF', # Upper Triassic\n", "          '#800080', # Middle Triassic\n", "          '#800080', # Lower Triassic\n", "          '#800080', # Lower Triassic\n", "         ]"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["fig = plt.figure(figsize=(15,4))\n", "plt.bar(df_wave.index,\n", "        df_wave['Wavelength'],\n", "        color=colors,\n", "        zorder=2,\n", "        width=0.5,\n", "       )\n", "plt.title('Example fold wavelengths for northern Kurdish stratigraphy\\nfrom Ramberg-Biot equation estimates',\n", "          fontsize=16\n", "         )\n", "plt.ylabel('Wavelength [m]', fontsize=14)\n", "plt.xticks(range(8), (wavelengths.keys()), rotation=90, fontsize=12)\n", "plt.grid(axis='y', zorder=1)\n", "plt.show()"]}, {"cell_type": "markdown", "metadata": {}, "source": ["### Geophysical equations\n", "\n", "If you're into geophysics, for many equations you won't even have to code them yourself as there is a ready-made repo called [bruges](https://github.com/agile-geoscience/bruges) which stands for:\n", "\n", "<br>\n", "<a href=\"https://agilescientific.com/bruges\" target=\"_blank\">\n", "    <img src=\"http://ageobot.elasticbeanstalk.com/bruges.png\" width=600 />\n", "</a>\n", "\n", "There are multiple modules inside Bruges:"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["import bruges\n", "\n", "dir(bruges)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["For example if you want to make a ricker wavelet, bruges has you covered:"]}, {"cell_type": "markdown", "metadata": {}, "source": ["<div style=\"padding:0.5em 1em 1em; background-color:#eeffee; border-radius:5px; border:1px solid #77aa77;\">\n", "\n", "<h2>EXERCISE</h2>\n", "\n", "<p><b>Use the bruges library to plot a ricker wavelet using matplotlib's `plt.plot()` function.</b></p>\n", "</div>"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["import ipywidgets as widgets\n", "\n", "@interact(duration=widgets.FloatSlider(value=0.256,min=0.04,max=0.512,step=0.004,\n", "                                      description='duration',\n", "                                      continuous_update=False,\n", "                                      readout_format='.3f'),\n", "          dt=widgets.FloatSlider(value=0.001,min=0.0001,max=0.008,step=0.0001,\n", "                                      description='dt',\n", "                                      continuous_update=False,\n", "                                      readout_format='.4f'),\n", "          frequency=widgets.IntSlider(value=25,min=1,max=75,step=1,\n", "                                      description='frequency',\n", "                                      continuous_update=False,\n", "                                      readout_format='d'),       \n", "          filled=widgets.Checkbox(value=True,description='fill wavelet',disabled=False)\n", "         )\n", "def plot_ricker(duration, dt, frequency, filled):\n", "    \"\"\"\n", "    Plot a filter:\n", "    Args:\n", "        function (function): one of ['ricker', 'gabor', 'sinc', 'cosine']\n", "        duration (float): The length in seconds of the wavelet.\n", "        dt (float): The sample interval in seconds.\n", "        frequency (ndarray): Dominant frequency of the wavelet in Hz.\n", "        fill (boolean): whether the filter plot is filled between 0 and wavelet.\n", "    Returns:\n", "        ndarray. {function} wavelet with centre frequency 'frequency' sampled on t.\n", "    \"\"\"\n", "    # call the wavelet function\n", "    w, t = bruges.filters.ricker(duration, dt, f=frequency, return_t=True)\n", "\n", "    # create the plot\n", "    fig, ax = plt.subplots(figsize=(15, 6), ncols=1) \n", "    ax.plot(t, w, color='black')\n", "    ax.grid()\n", "    ax.set_title(f'Ricker, frequency={frequency}, duration={duration}, dt={dt}')    \n", "    \n", "    # define fill_between() parameters\n", "    x_min = -duration / 2\n", "    x_max = duration / 2\n", "    x = np.arange(x_min, x_max, dt)\n", "    \n", "    # fill wavelet\n", "    if filled:\n", "        ax.fill_between(x, 0, w, where=w > 0, color='k')\n", "\n", "    # show the plot\n", "    plt.show()\n", "    \n", "    return\n"]}, {"cell_type": "markdown", "metadata": {}, "source": ["## Plotting a DEM and making cross sections\n", "In this section, we'll make a digital elevation map of Mount St-Helens after the [1980 eruption](https://en.wikipedia.org/wiki/1980_eruption_of_Mount_St._Helens) to see the crater left by the eruption, and we'll also plot some cross-sections through the mountain comparing the topography pre- and post-eruption.\n", "\n", "First we'll use [`numpy`](https://numpy.org/) to load some data we prepared earlier:"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["before = np.load('./data/st-helens_before.npy')\n", "after = np.load('./data/st-helens_after.npy')"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Then set up some variables:"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["difference = before - after\n", "fontsize = 14\n", "linewidth = 3\n", "ew_section1 = 190\n", "ns_section1 = 150"]}, {"cell_type": "markdown", "metadata": {}, "source": ["And finally make the [`figure` and `axes`](https://matplotlib.org/3.1.1/gallery/showcase/anatomy.html), first setting up the `figure` and then adding `axes` and specifying their positions to set up our plot just as we want it, we'll also add some `titles` and `labels`:"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["fig = plt.figure(figsize = (12,8))\n", "ax = fig.add_axes([0.1,0.1,0.8,0.8])\n", "ax.set_title('Mount St-Helens Digital Elevation Model map', fontsize=fontsize)\n", "ax.imshow(after)\n", "ax.axhline(ew_section1, color='red', lw=linewidth)\n", "ax.axvline(ns_section1, color='blue', lw=linewidth)\n", "ax.set_ylabel('Northing [index]', size=fontsize)\n", "ax.set_xlabel('Easting [index]', size=fontsize)\n", "ax.grid(alpha=0.5)\n", "\n", "ax2 = fig.add_axes([0.8,0.58,0.5,0.32])\n", "ax2.set_title('Mount St-Helens West-East cross-section', fontsize=fontsize)\n", "ax2.plot(before[ew_section1], c='red')\n", "ax2.plot(after[ew_section1], c='red')\n", "ax2.set_xlabel('Easting [index]', size=fontsize)\n", "ax2.set_ylabel('Elevation [m]')\n", "ax2.grid()\n", "\n", "ax3 = fig.add_axes([0.8,0.1,0.5,0.32])\n", "ax3.set_title('Mount St-Helens North-South cross-section', fontsize=fontsize)\n", "ax3.plot(before[:,ns_section1], c='blue')\n", "ax3.plot(after[:,ns_section1], c='blue')\n", "ax3.set_xlabel('Northing [index]', size=fontsize)\n", "ax3.set_ylabel('Elevation [m]')\n", "ax3.grid()"]}, {"cell_type": "markdown", "metadata": {}, "source": ["## Loading up well files and plotting curves\n", "\n", "First let's see what files we have with the inbuilt [`glob`](https://docs.python.org/3.7/library/glob.html#module-glob) module:"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "markdown", "metadata": {}, "source": ["### Inspecting a single well"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "markdown", "metadata": {}, "source": ["### Inspecting a single curve"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "markdown", "metadata": {}, "source": ["### Smoothing a curve"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "markdown", "metadata": {}, "source": ["### Plot all curves in one well"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "markdown", "metadata": {}, "source": ["### ...or in a group of wells"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "markdown", "metadata": {}, "source": ["### A Well has a `.df()` method that returns a pandas DataFrame\n", "\n", "We'll use the `.colums` attribute on the `DataFrame` to choose the logs we want:"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "markdown", "metadata": {}, "source": ["And then get some statistics for two curves using the [`.describe`](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.describe.html) method on the DataFrame:"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "markdown", "metadata": {}, "source": ["### Log QC\n", "\n", "The [`welly`](https://github.com/agile-geoscience/welly) library gives us access to some QC capabilities to rapidly inspect our logs - this might be particularly useful to prepare a Machine Learning task for example:"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["import welly.quality as qc\n", "from IPython.display import HTML\n", "\n", "tests = {\n", "    'each': [qc.no_gaps, qc.no_monotonic],\n", "    'GR': [qc.all_positive],\n", "    'RHOB': [qc.all_positive],\n", "}\n", "\n", "results = well_R39.qc_table_html(tests=tests)\n", "\n", "HTML(results)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["### Plot multiple logs of the same well\n", "\n", "Now we'll set up a log panel by using [`striplog`'s](https://github.com/agile-geoscience/striplog) [`Legend` module](https://striplog.readthedocs.io/en/latest/api/striplog.legend.html#module-striplog.legend). \n", "We select some logs we want to plot in different tracks, then define a legend for the display of each curve:"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["from striplog import Legend\n", "\n", "tracks = ['MD', 'CALI', 'GR', 'RHOB', 'NPOR', ['DT1R', 'DT2', 'DT2R']]\n", "\n", "legend = Legend.builtin('NSDOE')\n", "\n", "curve_legend_csv = \"\"\"colour,lw,ls,xlim,xscale,curve mnemonic\n", "black,0.5,-,,linear,CALI\n", "green,1.0,-,\"0,200\",linear,GR\n", "red,1.0,-,,linear,RHOB\n", "blue,1.0,-,,linear,NPOR\n", "red,0.4,-,,linear,DT1R\n", "#e6153b,0.4,-,,linear,DT2\n", "#c7324d,0.4,-,,linear,DT2R\n", "\"\"\"\n", "\n", "curve_legend = Legend.from_csv(text=curve_legend_csv)\n", "\n", "complete_legend = legend + curve_legend"]}, {"cell_type": "markdown", "metadata": {}, "source": ["And once this set up is done, we can make a plot for any given well using the `.plot` method and passing the `tracks` and `legend` variables to it:"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "markdown", "metadata": {}, "source": ["## Loading a seismic volume and inspecting slices\n", "\n", "For this task, we will use Equinor's [segyio](https://github.com/equinor/segyio) library which describes itself as follows:\n", "\n", "> Segyio is a small LGPL licensed C library for easy interaction with SEG Y formatted seismic data, with language bindings for Python and Matlab. Segyio is an attempt to create an easy-to-use, embeddable, community-oriented library for seismic applications. Features are added as they are needed; suggestions and contributions of all kinds are very welcome.\n", "\n", "It allows us to load 2D and 3D seismic, but of course with all the [usual caveats about the segy format](https://agilescientific.com/blog/2014/3/26/what-is-seg-y.html).\n", "\n", "So let's first import and inspect a curated seismic cube:\n", "\n", "If you don't have the file yet, **[get the large dataset from Agile's S3 bucket](https://s3.amazonaws.com/agilegeo/Penobscot_0-1000ms.sgy.zip)**. It's 140MB."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["import segyio\n", "\n", "with segyio.open('./data/Penobscot_0-1000ms.sgy') as s:\n", "    c = segyio.cube(s)"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "markdown", "metadata": {}, "source": ["<div style=\"padding:0.5em 1em 1em; background-color:#eeffee; border-radius:5px; border:1px solid #77aa77;\">\n", "\n", "<h2>EXERCISE</h2>\n", "\n", "<p><b>Plot a single section of the seismic, for example an inline section.</b></p>\n", "    <p>You'll have to make sure inlines and crosslines are <a href=\"https://numpy.org/doc/stable/reference/generated/numpy.transpose.html\">transposed</a> so that the seabed is horizontal.</p>\n", "</div>"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}, {"cell_type": "markdown", "metadata": {}, "source": ["## Scanning through a seismic volume interactively\n", "\n", "As we did before, we can use `interact` to explore this cube on the fly, as long as it fits in memory of course. We will define a `seismic loading` function, a `seismic plotting` function, and finally call this plotting function with `interact`:"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["def load_seismic(volume):\n", "    with segyio.open(volume) as s:\n", "        vol = segyio.cube(s)\n", "    return vol"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["def seismic_plotter(colormap, section, inline, xline, timeslice, volume):\n", "    \"\"\"Plot a given seismic ILine, XLine or Timeslice with a choice of colormaps\"\"\"\n", "    # load a seismic volume\n", "    vol = load_seismic(volume)\n", "    \n", "    # display options for different sections\n", "    sections = {\n", "        'inline': {'amp': vol[inline,:,:].T, 'line': inline, 'shrink_val': 0.6, \n", "                  'axhline_y': timeslice, 'axhline_c': 'b', \n", "                  'axvline_x': xline, 'axvline_c': 'g',\n", "                  'axspine_c': 'r',\n", "                  'aspect': 1\n", "                  },\n", "        'xline': {'amp': vol[:,xline,:].T, 'line': xline, 'shrink_val': 0.5, \n", "                  'axhline_y': timeslice, 'axhline_c': 'b', \n", "                  'axvline_x': inline, 'axvline_c': 'r',\n", "                  'axspine_c': 'g',\n", "                  'aspect': 1\n", "                 },\n", "        'timeslice': {'amp': vol[:,:,timeslice], 'line': timeslice, 'shrink_val': 0.95, \n", "                  'axhline_y': xline, 'axhline_c': 'g', \n", "                  'axvline_x': inline, 'axvline_c': 'r',\n", "                  'axspine_c': 'b',\n", "                  'aspect': 0.5\n", "                 },\n", "    }\n", "\n", "    # scale amplitudes\n", "    ma = np.percentile(vol, 98)\n", "    \n", "    # plot figure\n", "    fig, ax = plt.subplots(figsize=(18, 6), ncols=1)\n", "\n", "    # apply options\n", "    sec = sections[section]    \n", "    im = ax.imshow(sec['amp'], aspect=sec['aspect'], vmin=-ma, vmax=ma, cmap=colormap)\n", "    ax.set_title(f'Penobscot_0-1000ms {section} {sec[\"line\"]}', fontsize=18)\n", "    plt.colorbar(im, ax=ax, shrink=sec['shrink_val']).set_label('Amplitudes')\n", "     \n", "    # add projected lines\n", "    ax.axhline(y=sec['axhline_y'], linewidth=1.5, color=sec['axhline_c'])\n", "    ax.axvline(x=sec['axvline_x'], linewidth=1.5, color=sec['axvline_c'])\n", "    for axis in ['top','bottom','left','right']:\n", "        ax.spines[axis].set_linewidth(1.5)     \n", "        ax.spines[axis].set_color(sec['axspine_c'])\n", "    \n", "    return    "]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["_ = interact(seismic_plotter,\n", "             colormap=['Greys', 'Greys_r', 'seismic', 'seismic_r', 'viridis', 'plasma', 'inferno', 'magma'],\n", "             section=widgets.RadioButtons(options=['inline', 'xline', 'timeslice'],\n", "                                          value='inline',description='slicer',disabled=False),\n", "             inline=widgets.IntSlider(value=300,min=0,max=600,step=1,continuous_update=False,\n", "                                      description='<font color=\"red\">inline</>'),\n", "             xline=widgets.IntSlider(value=240,min=0,max=480,step=1,continuous_update=False,\n", "                                     description='<font color=\"green\">xline</>'),\n", "             timeslice=widgets.IntSlider(value=125,min=0,max=250,step=1,continuous_update=False,\n", "                                         description='<font color=\"blue\">timeslice</>'),\n", "             volume='./data/Penobscot_0-1000ms.sgy',\n", "            )"]}, {"cell_type": "markdown", "metadata": {}, "source": ["<hr />\n", "<img src=\"https://avatars1.githubusercontent.com/u/1692321?v=3&s=200\" style=\"float:center\" width=\"40px\" />\n", "<p><center>\u00a9 2020 <a href=\"http://www.agilegeoscience.com/\">Agile Geoscience</a> \u2014 <a href=\"https://creativecommons.org/licenses/by/4.0/\">CC-BY</a></center></p>"]}], "metadata": {"kernelspec": {"display_name": "transform2020", "language": "python", "name": "transform2020"}, "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": 4}