t21-segysak.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "forty-easter",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<img src=\"images/swung_logo_vector.png\" alt=\"swung\" style=\"width: 40%;\"/>\n",
    "\n",
    "# T21 Segysak Tutorial - Tony Hallam, April 2021\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "charged-scratch",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "Welcome!\n",
    "\n",
    "**Admin:**\n",
    "  - Get the tutorial and data on Github (https://github.com/trhallam/segysak-t21-tutorial)\n",
    "  - Run the tutorial on binder.\n",
    "  - Limited interaction during the video but talk to the `segysak` experts in GatherTown and on the Tutorial Slack Channel"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "affiliated-product",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Introduction\n",
    "\n",
    "## What is `segysak`?\n",
    "\n",
    "## Tour\n",
    "\n",
    "## `segysak` vs X\n",
    "\n",
    "## SEG-Y Files"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "outdoor-mixer",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## `segysak`\n",
    "\n",
    "Make SEG-Y data easily accessible and creatable from Python"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "specified-recorder",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "Simply put, `segysak` has grown as a set of tools to make SEG-Y data easily accessible and createable from Python.\n",
    "It leverages a number of existing libraries but brings them together to try and improve the user experience, and\n",
    "to remove as much boiler plate code as possible when dealing with SEG-Y.\n",
    "\n",
    "The project started about a year ago at Transform 2020. Most of the work was done during that hackathon but it has\n",
    "continued to develop since then with gradual improvements, bug-fixes and user support.\n",
    "\n",
    "Although I'm the project owner and one of the primary users of `segysak` (I use it it in a lot of my PhD projects).\n",
    "It is open for the subsurface community to not only utilise, but to contribute to grow to meet peoples needs.\n",
    "\n",
    "I'd strongly encourage anyone with ideas and/or enthusiasm for changes or additions to get in touch so we can improve `segysak` for everyone."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "infrared-titanium",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Tour\n",
    "\n",
    " - Github (source code, issues, contributions) - https://github.com/trhallam/segysak\n",
    " - Documentation (help, examples, API) - https://segysak.readthedocs.io/en/latest/\n",
    " - Slack (help, discussion, ideas, contributions) - https://swung.slack.com/messages/segysak/"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "intense-requirement",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "Everything you need to know about `segysak` is available online. There is the Github repository where we manage the source code for the library and distribute the packages for installation via pip. \n",
    "\n",
    "There is also an issue tracker where you can raise bugs/problems or submit ideas or suggestions. It's also a good place to look for things that need doing if you want to help out.\n",
    "\n",
    "We then have the documentation on readthedocs. Here you will find more detailed help, examples (which are avaialable as Notebooks) and the API (of function and member descriptions). This is a really useful place to come if you are stuck, or\n",
    "need more detail because we cover a lot of the basics in the documentation. Indeed this workshop is heavily influenced\n",
    "by the first few example notebooks you can find here.\n",
    "\n",
    "Finally, we have the Slack forum hosted on swung.slack. This space is always open for people to ask questions or get help, even drop by just for a bit of discussion."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "norman-minute",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## `segysak` versus X"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "pursuant-husband",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "A lot of the time I get asked about segysak versus X in the Python world, where does it fit in?\n",
    "The reality is, segysak doesn't so much compete with any part of the scientific stack but tries to form bridges over\n",
    "the common space we often have to traverse. For example."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "amazing-compensation",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "### `segyio`\n",
    "\n",
    " - `segysak` relies on `segyio` but abstracts a lot of the low level detail"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "previous-unemployment",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "segysak couldn't exist without segyio - segyio does all the hard work of interacting with the actual SEG-Y and segysak tries to make segyio a bit more accessible by providing a direct link between it and easy to use libraries like xarray."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "quantitative-concrete",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "### `xarray`\n",
    "\n",
    " - `segysak` extends `xarray` to make it easier to deal with SEG-Y files"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dramatic-nursery",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "Things like loading and writing of files are more automated. Trys to take care of tracking things like headers, and attributes for you.\n",
    "\n",
    "Also includes extensions for common seismic related tasks."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tamil-pontiac",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## SEG-Y Files\n",
    "\n",
    "File format defined by the SEG Organisation for storing seismic trace data.\n",
    "\n",
    "Heavily geared toward limited size magnetic reel tapes."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abstract-theory",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "SEG-Y data is a legacy format developing during a time when storage space was significantly more limited than the present. Specifically it was designed to allow the recording of shot records onto magnetic tapes. Whilst it has generally served the industry well, it is not particularly performant and there are a lot of inconsistencies between what information is included with a file and where is stored.\n",
    "\n",
    "The most consistent aspects of a SEG-Y file are consistent though, and that is size and layout of the repeated header and data sections of the file. \n",
    "\n",
    "With `segysak` there is functionality to allow you to read some or all of the data and to inspect or extract information from the various parts of the headers."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "headed-radius",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Basic Format (SEGY-Rev2):**\n",
    "\n",
    "<img src=\"images/segy_layout.png\" alt=\"swung\" style=\"width: 100%;\"/>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "subject-large",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Installation and Tutorial Setup\n",
    "\n",
    "`pip install segysak`\n",
    "\n",
    "Demo Data\n",
    "\n",
    "Opening the Tutorial Notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "personal-market",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "`segysak` is currently only available from source or via `pip`.\n",
    "\n",
    "To run this tutorial locally, you'll need to clone the entire tutorial repository which also contains the demo data. It is then just a matter of firing up Jupyter and opening the tutorial notebook.\n",
    "\n",
    "We can check to make sure everything is running by importing the `segy_loader` function from `segysak.segy`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "sonic-teaching",
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from segysak.segy import segy_loader\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "private-grace",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Basic Imports and Test Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "southwest-offering",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import pathlib\n",
    "import platform\n",
    "from IPython.display import display\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "alpha-schema",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# specify the example file and check we have the example data\n",
    "\n",
    "segy_file = pathlib.Path(\"data/volve10r12-full-twt-sub3d.sgy\")\n",
    "print(\"SEG-Y exists:\", segy_file.exists())\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "professional-response",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Basic Usage\n",
    "\n",
    "### Inspecting SEG-Y files\n",
    "\n",
    "### Loading SEG-Y files\n",
    "\n",
    "### `xarray.Dataset` basics\n",
    "\n",
    "### NetCDF Files\n",
    "\n",
    "### Editing and Saving"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "characteristic-andrews",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Inspecting SEG-Y files\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "composite-forge",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "There are a number of utility functions in `segysak` designed to help you explore and understand the data in SEG-Y files without needing to load the entire SEG-Y file in.\n",
    "\n",
    "For example we can look at the text header of the file by using the function `get_segy_texthead`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "confused-youth",
   "metadata": {},
   "outputs": [],
   "source": [
    "from segysak.segy import segy_header_scan, segy_header_scrape, get_segy_texthead"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "angry-testing",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# examine the text header\n",
    "get_segy_texthead(segy_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "selective-string",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Inspecting SEG-Y files - trace header scan\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "verified-latitude",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "The `segy_header_scan` function allows us to get a `pandas.DataFrame` containing information about the contents of the headers for the first few traces of a SEG-Y file. This saves us having to scan the whole file and often contains enough information for us to then load the file properly.\n",
    "\n",
    "The DataFrame index names are the same as what `segyio` uses but `segysak` tables them and gives you a few bits of information. Importantly for properly loading SEG-Y data you will need to note the byte locations of header information you want. \n",
    "\n",
    "You can increase the number of traces you want to scan by setting the `max_traces_scan` keyword in the function - by default it is 1000 traces.\n",
    "\n",
    "Try using the context manager to display more rows.\n",
    "```python\n",
    "with pd.option_context(\"display.max_rows\", 100):\n",
    "    display(df)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "expressed-litigation",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# scan the headers to check\n",
    "scan = segy_header_scan(segy_file, max_traces_scan=2000)\n",
    "scan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "numeric-canyon",
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "with pd.option_context(\"display.max_rows\", 100):\n",
    "    display(scan)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "intimate-updating",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Inspecting SEG-Y files - trace header scrape\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "oriented-ministry",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "We can also extract all of the trace header information using `segy_header_scrape`. This function creates a complete copy of the traces headers as a `pandas.DataFrame`. On larger files it can take a little bit of time to scan all the traces. You can see even on this small volume there are 12,322 traces.\n",
    "\n",
    "The column names again are the same as in `segyio`.\n",
    "\n",
    "We can do a quick check of the headers here by creating some simple plots. If you're used to loading SEG-Y in commercial software they usally offer you something showing trace number vs value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tired-metadata",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "trace_headers = segy_header_scrape(segy_file)\n",
    "trace_headers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "written-fault",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(nrows=2, ncols=2, figsize=(20, 10), sharex=True)\n",
    "\n",
    "for ax, prop in zip(axs.ravel(), [\"CDP_X\", \"CDP_Y\", \"INLINE_3D\", \"CROSSLINE_3D\"]):\n",
    "    ax.plot(trace_headers[prop])\n",
    "    ax.set_title(prop, fontdict={\"fontsize\":18})\n",
    "\n",
    "for ax in axs[1, :]:\n",
    "    ax.set_xlabel(\"Trace\", fontdict={\"fontsize\":18})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "printable-desperate",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Loading SEG-Y files - complete\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cardiac-moldova",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "Loading of SEG-Y data is also pretty straight forward, and in `segysak` there are a few different ways to go about it.\n",
    "\n",
    "The most straightforward way is to use the multi-purpose function `segy_loader`. This function loads, 2D, 3D and gathers. There are a lot of options to customise the loader for most situations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tough-machinery",
   "metadata": {},
   "outputs": [],
   "source": [
    "from segysak.segy import segy_loader\n",
    "help(segy_loader)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "alike-traffic",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# loading with default byte locations\n",
    "seisnc_vol = segy_loader(segy_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "humanitarian-burlington",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "`segy_loader` returns an `xarray.Dataset` which has dimensions appropriate for the type of seismic loaded. In this case we have a 3D volume so our dimensions are `iline`, `xline`, and because we used the default option of `TWT` the vertical dimension is `twt`, but if you have a depth volume you can specify this in the `segy_loader` keyword arguments using `vert_domain=\"DEPTH\"`.\n",
    "\n",
    "The actual seismic data/volume is contained within the `data` variable, and we see this has the full dimensions of the cube. Whilst the `cdp_x` and `cdp_y` values from the trace headers don't have the vertical dimension.\n",
    "\n",
    "There is also a number of attributes which are created by the loader and provide information about the loaded file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "random-spine",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(seisnc_vol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "concerned-genealogy",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# lets quickly checkout our data - we'll talk about xarray basics in the next section\n",
    "_ = seisnc_vol.sel(iline=10100).data.T.plot(yincrease=False, figsize=(20, 10), vmax=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "narrative-horror",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "Here is an example where we explicitly set the key header byte locations of `iline`, `xline`, `cdpx` & `cdpy`. If other values from the header are needed, the byte locations can be set using the `extra_byte_fields` keyword argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "saving-arrangement",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# specifying byte locations for key cube geometry\n",
    "_ = segy_loader(\n",
    "    segy_file,\n",
    "    iline=189, xline=193, cdpx=181, cdpy=185,\n",
    "    vert_domain=\"DEPTH\",\n",
    "    #extra_byte_fields=[117]\n",
    "    extra_byte_fields={\"my_name\":117}\n",
    ")\n",
    "print(_)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "prostate-traveler",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Loading SEG-Y files - with filtering\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "filled-massachusetts",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "It is also possible to filter the data being loading in using the cropping keywords. These can be used to restrict the amount of data being loaded in. In this example we load just a single in-line `10100`. \n",
    "\n",
    "Unfortunately to filter the data we still have to scan all the headers but the loaded volume is now smaller."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "republican-abortion",
   "metadata": {},
   "outputs": [],
   "source": [
    "seisnc_vol_iline_10100 = segy_loader(segy_file, ix_crop=(10100, 10100, 2000, 3000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "surgical-shaft",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(seisnc_vol_iline_10100)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "satellite-shock",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "Sometimes SEG-Y are really big and scanning the headers repeatedly can be awkward. \n",
    "Recall that we previous scanned the whole trace headers into a DataFrame. \n",
    "`segysak` allows us to perform filtering on the DataFrame, and use the filtered version for loading. \n",
    "This provides an enormous amount of flexibility on what data is loaded and means the headers only need to be scanned once on multiple loading events.\n",
    "\n",
    "The trace header dataframe should be passed to the `head_df` keyword argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "prescription-power",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "seisnc_vol_iline_block = segy_loader(\n",
    "    segy_file,\n",
    "    head_df=trace_headers[trace_headers[\"INLINE_3D\"] <= 10100].copy()\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "coupled-chuck",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(seisnc_vol_iline_block)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "editorial-generation",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Other useful loading functions in `segysak`.\n",
    "  - `segy.segy_converter` - Streams data from SEG-Y to NetCDF4 on disk\n",
    "  - `segy.segy_freeloader` - Support for higher dimensional data (development branch).\n",
    "  - `openzgy.zgy_loader` - Experimental support for ZGY based upon open ZGY."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "editorial-appointment",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## `xarray.Dataset` basics\n",
    "\n",
    "Based upon the NetCDF file format for multi-variable, n-dimensional data.\n",
    "\n",
    "<center> <img src=\"images/seisnc-diagram.png\" alt='seisnc' width=\"50%\" /> </center>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "italian-offering",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "All of the intricacies of `xarray` are a bit beyond this tutorial but we'll try to quickly cover here some of the most useful ones with `segysak` seismic Datasets.\n",
    "The different parts of the dataset can be accessed through properties of the class.\n",
    "  - dimensions : `dim`\n",
    "  - coordinates : `coords`\n",
    "  - variables : `variables`\n",
    "  - attributes : `attrs`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "special-possession",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# dataset anatomy - dimensions, coordinates, variables, attributes - DataArray vs Dataset\n",
    "seisnc_vol"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "resistant-inventory",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Sub-selections with `xarray` are not made by using regular indexing like `numpy` for example because `xarray` does not guarantee the order of dimensions. The key methods for selection are `sel`, and `isel` which allow labelled selection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "hindu-hierarchy",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# data selection - sel, isel"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "moral-exchange",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The values of any dimension, coordinates, and variables can be returned using the property."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "continuous-philadelphia",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# data as numpy array\n",
    "seisnc_vol.iline.values"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "modified-manner",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "There are also a few useful methods to know about including `plot`, `interp`, `mean`, `max`, `min` and so on.\n",
    "\n",
    "Two really important ones are `transpose` which lets you do labeled transposing of variables and `broadcast_like` which allows you to transform one variable or dataset to match another. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "consolidated-difference",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# xarray methods (plot, interp, mean, min, max, transpose, broadcast_like, etc...)\n",
    "# seisnc_vol.cdp_x.transpose(\"xline\", \"iline\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "gorgeous-fusion",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The last important thing to learn with `xarray` is how to assign data into variables. Because `DataSets` are multi-dimensional we also have to give `xarray` information about the dimensions of the data.\n",
    "\n",
    "```python\n",
    "print(seisnc_vol.dims)\n",
    "seisnc_vol[\"zeros\"] = ((\"iline\", \"xline\", \"twt\"), np.zeros((61, 202, 850)))\n",
    "seisnc_vol\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "least-green",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# xarray variable assignment"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "hawaiian-topic",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## `xarray` FAQ\n",
    "\n",
    " - Why don't we make the global coordinates the dimensions?\n",
    " - How do I save/persist my changes."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "afraid-basis",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "Notes: \n",
    "Global coordinates are not orthogonal because the seismic grid rarely lines up with Grid North.\n",
    "Persisting changes either means saving back to SEG-Y or using the NetCDF File Format. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "marked-stereo",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# NetCDF File Format\n",
    "\n",
    "Common in climate science, binary, fast and lazy loading\n",
    "\n",
    "(basically `xarray.Datataset` on disk)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "faced-helen",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## NetCDF\n",
    "\n",
    "Why use “another” file format for seismic?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fluid-peoples",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "Generally it just makes working with seismic in Python easier. It will save you time if you are reading volumes repeatedly or can't store everything you need in memory.\n",
    "\n",
    "NetCDF was the logical choice because it is at the core of `xarray` but `xarray` supports other data models such as zarr which we are investigating.\n",
    "There is also beta support within `segysak` for the OpenZGY format with instructions on Github about how to set that up."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "announced-worst",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    " - Faster than SEG-Y for most use cases."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "positive-northwest",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    " - Widely supported within the Python scientific stack (xarray, dask)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "offshore-fourth",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    " - Commonly supported in other languages."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "genuine-lying",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Saving the data to netcdf requires the use of the seisio accessor due to limitations on the types of attributes that can be\n",
    "saved using the xarray method.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "arabic-mount",
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# output the data to netcdf\n",
    "seisnc_vol.seisio.to_netcdf(\"data/test.seisnc\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dimensional-saying",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "if platform.system() == \"Windows\":\n",
    "    !dir data\\.\n",
    "else: ## linux\n",
    "    !ls data/."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "referenced-globe",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We can check to see if the seisnc NetCDF file was created ok by reimporting it. To get it back into the same form that `segysak` uses we can use the `open_seisnc` function. `open_seisnc` is a thin wrapper around the `xarray.open_dataset` method that includes some special handling for segysak attributes and ensures that the dataset is opened with the `.seis` extension for xarray which we will talk about soon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "naughty-engineering",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from segysak import open_seisnc\n",
    "open_seisnc(\"data/test.seisnc\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "clear-dallas",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Saving to SEG-Y"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "native-natural",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "Generally if you have loaded a SEG-Y file and edited it, you can then save that file back to a new SEG-Y in one line.\n",
    "Currently SEGY-SAK doesn't support editing SEG-Y in place but it is something that might come in the future if the demand is\n",
    "there (or you could help develop this for us). \n",
    "\n",
    "There are a few attributes that `segysak` needs to write your new SEG-Y file.\n",
    "  - `coord_scalar` (int)\n",
    "  - `sample_rate` (float)\n",
    "  - `source_file` (str)\n",
    "  \n",
    "It also needs dimensions from one of the dimension sets. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "equal-investment",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from segysak.segy import segy_writer\n",
    "help(segy_writer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ancient-ghost",
   "metadata": {},
   "outputs": [],
   "source": [
    "test = seisnc_vol.copy()\n",
    "test.attrs = {\"coord_scalar\":-100, \"sample_rate\":4.0, \"source_file\":\"\"}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "boring-operation",
   "metadata": {},
   "source": [
    "If byte locations need to be changed for specific software the `trace_header_map` keyword is available. Any variable keys in the dataset can be assigned to a trace header byte location.\n",
    "\n",
    "In this example we set the `iline` to go to byte location 21."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "moral-dress",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# export in memory dataset to segy\n",
    "segy_writer(test, \"data/test.segy\", trace_header_map={\"iline\":21})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "refined-yeast",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "If we create variables that cover the trace header dimensions (iline/xline) then these can also be included in the output to SEG-Y by specifying the variable key and the byte location where the variable should be placed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "saving-korea",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# write other variables\n",
    "seisnc_vol[\"xy\"] = seisnc_vol[\"cdp_x\"] * seisnc_vol[\"cdp_y\"] / 1E10\n",
    "seisnc_vol\n",
    "# segy_writer(seisnc_vol, \"data/test.segy\", trace_header_map={\"xy\":13})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "starting-franchise",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Lets read the headers of that SEG-Y and see what we got. The output has been converted to int so our floating point values are gone. This is a limitation of the SEG-Y format. Any floating point numbers must be scaled to int and back again on loading. This is done automatically for coordinates but all other values must be handled manually."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "anticipated-patrol",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "segy_header_scrape(\"data/test.segy\").T"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "industrial-water",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 10 Minute Break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "nuclear-annual",
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "from tqdm.auto import tqdm\n",
    "\n",
    "with tqdm(desc=\"Break Timer\", total=10*60, bar_format=\"{l_bar}{bar} {elapsed_s:.0f}/{total} seconds\") as pbar:\n",
    "    start = time.time()\n",
    "    now = time.time()\n",
    "    prev_now = now\n",
    "    while (now - start) < 10*60:\n",
    "        pbar.update(now - prev_now)\n",
    "        time.sleep(1)\n",
    "        prev_now = now\n",
    "        now = time.time()\n",
    "    pbar.update(time.time() - prev_now)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dress-burner",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Horizon extraction\n",
    "\n",
    " - Load a horizon and add it to a cube\n",
    " - Plotting maps\n",
    " - Plotting horizons on vertical slices\n",
    " - Extracting seismic ampltidues along a horizon"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "static-federal",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Load some seismic horizon data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "chinese-bibliography",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Lets start by specifying the path to some seismic horizon data and checking it is available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "equal-singer",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "top_hugin_path = pathlib.Path(\"data/hor_twt_hugin_fm_top.dat\")\n",
    "print(\"File\", top_hugin_path, \"exists?\", top_hugin_path.exists())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "supreme-holiday",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "If we quickly look at the first few lines of the file we can see it is a space delimited file with three columns.\n",
    "UTM X, UTM Y and TWT. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "metric-mileage",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# check the file layout\n",
    "with open(top_hugin_path) as f:\n",
    "    lines = [next(f) for i in range(5)]\n",
    "print(*lines)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "copyrighted-support",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "It is then quite straightforward to load the file in using `pandas.read_csv`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "starting-astronomy",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# is a csv file\n",
    "top_hugin_df = pd.read_csv(top_hugin_path, names=[\"cdp_x\",\"cdp_y\",\"twt_hugin\"], sep=' ')\n",
    "top_hugin_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "genetic-louisiana",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "When the horizon is in this format it might not map directly to the seismic trace locations. To simplify the process of interpolating the horizon to the seismic trace locations, `segysak` has a `surface_from_points` method in the `.seis` accessor. By default, this method will try to interpolate using `cdp_x` and `cdp_y` from the dataset but these options can be changed.\n",
    "\n",
    "A new dataset is returned with the same dimensions as the seismic volume but now with the horizon data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "hearing-olympus",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "top_hugin_ds = seisnc_vol.seis.surface_from_points(top_hugin_df, 'twt_hugin', right=('cdp_x', 'cdp_y'))\n",
    "print(top_hugin_ds)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "refined-launch",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Plotting Maps"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "municipal-mailman",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "For plotting we can use the built in plot command that comes with `xarray` datasets. This is wrapper around a call to `matplotlib` and allows us to get quickly formatted plots of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "seeing-collaboration",
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "top_hugin_ds.twt_hugin.plot(cmap='hsv')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aware-puppy",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "This plotting is done on the local iline/xline grid which defines the cube though, and often we want to see things in a UTM X and  Y type context.\n",
    "\n",
    "There are a couple of ways to go about this.\n",
    "  - One way is use explicit plotting based upon the x and y coordinates in the dataframe.\n",
    "  - Another is to use a transform argument for `matplotlib`'s plotting commands."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "intense-ocean",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "axs = plt.subplot()\n",
    "mesh = axs.pcolormesh(\n",
    "    top_hugin_ds.cdp_x.values,\n",
    "    top_hugin_ds.cdp_y.values,\n",
    "    top_hugin_ds.twt_hugin.values,\n",
    "    shading=\"auto\"\n",
    ")\n",
    "axs.set_aspect(1)\n",
    "_ = plt.colorbar(mesh, orientation=\"horizontal\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "solar-justice",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Using the transform can be useful when we want to plot objects using iline/xline notation but in x-y coordinate space, like an inline location for example (10100). We can also use the inverted form of the transform to convert x and y coordinates to iline and xline. `segysak` also has the `.seis.xysel` method to extract seismic based upon x and y trace locations. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "liquid-bruce",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "tform = seisnc_vol.seis.get_affine_transform()\n",
    "\n",
    "axs = plt.subplot()\n",
    "mesh = axs.pcolormesh(\n",
    "    top_hugin_ds.iline,\n",
    "    top_hugin_ds.xline,\n",
    "    top_hugin_ds.twt_hugin.T,\n",
    "    shading=\"auto\",\n",
    "    transform=tform + axs.transData\n",
    ")\n",
    "axs.set_aspect(1)\n",
    "_ = axs.plot([10100, 10100], [2200, 2300], transform=tform + axs.transData, color=\"w\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cheap-pulse",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Plotting Horizons on vertical section views"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "expanded-carol",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The horizon data can also be assigned back to the original seismic dataset. This can be useful for doing simultaneous sub-selection of the two variables at once. The dimensions of the data don't need to be specified in this case because the `top_hugin_ds` dataset already has the dimensions assigned. `xarray` will automagically line up those named dimensions with common named dimensions in `seisnc_vol`. Where the dimensions don't line up or new dimensions are included they will be extended or added to `seisnc_vol`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "english-exclusive",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# assign horizon back to seismic\n",
    "seisnc_vol[\"hugin\"] = top_hugin_ds.twt_hugin\n",
    "print(seisnc_vol)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "surprising-commander",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Combining the horizon into the seismic dataset can make many subsequent tasks much simpler. For example, here we sub-select a single inline once, and then use the reference in two subsequent plotting calls, once for the seismic, and then again for the horizon. Note that the state of `iline_subsel` as a view or copy depends on a lot of things that are a little hard to predict with `xarary`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "round-christianity",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# plotting\n",
    "iline_subsel = seisnc_vol.sel(iline=10100, twt=range(2402, 2900, 4), method='nearest')\n",
    "fig, axs = plt.subplots(figsize=(20, 5))\n",
    "iline_subsel.data.T.plot(ax=axs, yincrease=False)\n",
    "_ = axs.plot(iline_subsel.xline, iline_subsel.hugin, 'k')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "opposed-basketball",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Seismic amplitude maps"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "respected-cameroon",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Extracting the intersection of a horizon with a seismic volume is realy simple in `xarray`. It is literally one line. In this case, `xarray` understands the `iline` and `xline` relationship between the input DataArray `seisnc_vol.hugin` and the seismic volume. When we pass it via the `interp` method, `xarray` performs a linear interpolation to find the intersection point returning a new amplitude DataArray."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "patent-depression",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "amp = seisnc_vol.data.interp({\"twt\": seisnc_vol.hugin}, method='linear')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "miniature-huntington",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "axs = plt.subplot()\n",
    "mesh = axs.pcolormesh(amp.iline, amp.xline, amp.T, transform=tform + axs.transData, shading=\"auto\", cmap=\"bwr_r\", vmin=-6, vmax=6)\n",
    "ctr = axs.contour(top_hugin_ds.cdp_x, top_hugin_ds.cdp_y, top_hugin_ds.twt_hugin, colors='w')\n",
    "axs.set_aspect(1)\n",
    "plt.colorbar(mesh)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "imported-remark",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Mapping functions over blocks\n",
    "\n",
    " - Learn how to use Xarray to map functions on blocks of data, such as trace maths."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "applied-avatar",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "  - Horizon Flattening\n",
    "<img src=\"images/hflat.png\" alt=\"hflat\" style=\"width: 60%;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "antique-surprise",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "To do this we will introduce the `groupby` method for a dataset but to make it useful we need to create a trace identifier. `groupby` actually uses `pandas` in the backend and is the same. It will create a group of datasets based upon the key or keys that can be submitted to a function to perform operations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "about-argument",
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "for grp, subds in seisnc_vol.groupby(seisnc_vol.iline):\n",
    "    print(grp)\n",
    "    print(subds)\n",
    "    break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "expired-phenomenon",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# create a trace identifier\n",
    "seisnc_vol[\"trace\"] = ((\"iline\", \"xline\"), np.arange(61*202, dtype=int).reshape(61, 202))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "under-operation",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "for grp, subds in seisnc_vol.groupby(\"trace\"):\n",
    "    print(grp)\n",
    "    print(subds)\n",
    "    break"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "informal-former",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The next step is writing a function that takes advantage of the data as it is made available from `groupby`.\n",
    "To flatten the group be need to shift the time axis for each trace so that the horizon occurs at a constant time. The simplest way to do this is just to make that constant zero, so we subtract the horizon value form the time axis. \n",
    "\n",
    "Xarray also requires that our output cube be regularly sampled, so arbitraty shifts for each trace need to be resampled to a regular grid. That can be done using the `interp` function. `interp` applies a chosen interpolation methods (in this case linear interpolation) to resample the data against a new twt axis which we will call `twt_out`. Then the resampled trace is returned."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "statistical-ireland",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def hflat(ds, hor_var, twt_out):\n",
    "    trace_out = ds.copy()\n",
    "    trace_out[\"twt\"] = ds.twt - np.squeeze(ds[hor_var].values)\n",
    "    return trace_out.data.interp(twt=twt_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "supported-joshua",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We are also going to need to specify what the output time dimension should be (`twt_out` in our function). We know that subtracting the horizon TWT from the TWT grid will result in a new full TWT grid that goes from -hor_max to twt_max-hor_min of the input horizon. So we create a flattened TWT range with a sample interval of 1ms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "engaged-remark",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "flat_twt = np.arange(-seisnc_vol.hugin.max(), seisnc_vol.twt.max()-seisnc_vol.hugin.min(), 1, dtype=int)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "molecular-harvey",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "To recombine all the data back into a single `dataset` we can tag the `map` function onto the end of `groupby`. Doing this automatically applied what pandas and xarray call `split-apply-combine` logic, and is really handy.\n",
    "\n",
    "Here I'm just going to apply the process to a single inline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "about-district",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# applying groupby().map()\n",
    "tg_gby = seisnc_vol.sel(iline=10100).groupby(\"trace\").map(hflat, args=(\"hugin\", flat_twt))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "affecting-detroit",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Lets plot up our normal and flattened volumes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "traditional-gates",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# plotting results\n",
    "fig, axs = plt.subplots(ncols=2, figsize=(20, 5))\n",
    "seisnc_vol \\\n",
    "    .sel(iline=10100, twt=range(2002, 2900, 4), method='nearest') \\\n",
    "    .data.T.plot(ax=axs[0], yincrease=False)\n",
    "axs[0].plot(seisnc_vol.sel(iline=10100).xline, seisnc_vol.sel(iline=10100).hugin, 'k')\n",
    "tg_gby \\\n",
    "    .sel(twt=range(-300, 300, 4), method='nearest') \\\n",
    "    .T.plot(ax=axs[1], yincrease=False)\n",
    "axs[1].hlines(0, 0, 10000, \"k\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "static-bearing",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Groupby is generally fine for small volumes, but when you start to scale up your datasize it can run into issues. Xarray links nicely to the `dask` distributed processing library and it can even lazyily load data from disk that won't fit into memory. There is `dask` tutorial on the SEGY-SAK RTD website but today I'm just going to demonstrate how we can use the experimental `map_blocks` function to to achieve the same outcome as `groupby` but in a `dask` friendly way.\n",
    "\n",
    "Xarray has a chunking feature which allows you to break your Dataset down into smaller operational blocks. In this case we want the blocks to match the size of our operations which is just 1 trace. \n",
    "\n",
    "Because we won't necesarilly be doing the whole operation at once with `map_blocks` we also need to tell xarray what we think the output of our function will look like. This is done using a template. In our case the output is the same in every way, except the time dimension is resampled. Here we can just resample that dimension using `interp` to get the right output template shape.\n",
    "\n",
    "Note that the template is not actually calculated, just a place holder is created. This is because when we chunk xarray Datasets or DataArrays, every operation is delayed until the last possible moment. This is good for memory management and really simplifies the whole process. In this instance if you wanted the template computed you could either call `template.compute()` or access the numpy `values` array of the data variable.\n",
    "\n",
    "Then, similar to `groupby` we pass `map_blocks` the flattening function, the extra arguments, and now also the template."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "appropriate-software",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "seisnc_vol_chkd = seisnc_vol.chunk({\"iline\":1, \"xline\":1})\n",
    "template = seisnc_vol_chkd.interp(twt=flat_twt)\n",
    "tg_mb = seisnc_vol_chkd.map_blocks(hflat, args=(\"hugin\", flat_twt), template=template.data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "exact-cardiff",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "`tg_mb` is also a delayed object here. There is a task registered for each trace, but they will only be calculated when we ask for the data. Such as when we perform plotting. To get the full volume, we have to call `compute` again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bored-identity",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(ncols=2, figsize=(20, 5))\n",
    "seisnc_vol.sel(iline=10100, twt=range(2002, 2900, 4), method='nearest').data.T.plot(ax=axs[0], yincrease=False)\n",
    "tg_mb.sel(iline=10100, twt=range(-300, 300, 4), method='nearest').T.plot(ax=axs[1], yincrease=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "conventional-afternoon",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## `map_blocks` final thoughts\n",
    "\n",
    "  - Not super mature yet but very useful.\n",
    "  - Slightly different way of thinking to Python's normal instant run/result.\n",
    "  - Other useful delayed functions are: `to_netcdf`, `rolling`, `interp`, but most xarray operations can be delayed."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "incorporate-effect",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Vectorization of Seismic\n",
    "\n",
    " - I want to do machine learning and I need to tabularize my seismic and headers.\n",
    " - Now I need to send my results back to SEG-Y.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "competitive-blues",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Converting an `xarray.Dataset` to a `pandas.DataFrame` is really simple due to the close ties between the two packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cognitive-lawrence",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# creating a table from seismic\n",
    "seisnc_vol_df = seisnc_vol.isel(iline=10).to_dataframe()\n",
    "print(seisnc_vol_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "together-solution",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The Dataframe will have what `pandas` calls a multi-index, so to remove it we just need to reset the index. Note this will have a big impact upon your memory footprint."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "hispanic-samuel",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "seisnc_reindex = seisnc_vol_df.reset_index()\n",
    "print(seisnc_reindex)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "superb-property",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(seisnc_vol_df.info())\n",
    "print(seisnc_reindex.info())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "yellow-custody",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "When operations have been completed in the tabular format we can return to `xarray`. First the multi-index must be restored to get the coordinates right."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "devoted-tobago",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "seisnc_df_multi = seisnc_reindex.set_index([\"iline\", \"xline\", \"twt\"])\n",
    "print(seisnc_df_multi)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "spatial-tiffany",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "And then we just use the `to_xarray` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "amber-string",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "seisnc_xr = seisnc_df_multi.to_xarray()\n",
    "print(seisnc_xr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "micro-grant",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "The process isn't perfect though, we can see that 'cdp_x' and 'cdp_y' have come back as 3d cubes. And we will need to reset all the seisnc attributes are missing before we can export the data to SEG-Y using `segy_writer`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "animated-cameroon",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "seisnc_xr.attrs = seisnc_vol.attrs\n",
    "display(seisnc_xr.attrs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "banner-parameter",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "seisnc_xr[\"cdp_x\"] = seisnc_xr[\"cdp_x\"].mean(dim=[\"twt\"])\n",
    "seisnc_xr[\"cdp_y\"] = seisnc_xr[\"cdp_y\"].mean(dim=[\"twt\"])\n",
    "seisnc_xr = seisnc_xr.set_coords([\"cdp_x\", \"cdp_y\"])\n",
    "print(seisnc_xr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "controversial-possession",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Questions - Slack time because we will run over.\n",
    "\n",
    " - Fall backs to chat about memory management, dask, other file formats such as ZGY and Zarr\n",
    " - Demo of CLI for quick looks at headers or EBCIDC\n",
    " - Contribution Opportunities / Community led development\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "major-assumption",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (segysak)",
   "language": "python",
   "name": "segysak"
  },
  "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.6"
  },
  "rise": {
   "backimage": "segysak_background.png"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}