{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "db000e60009e44a4",
   "metadata": {},
   "source": [
    "# Read Data from the Cloud\n",
    "\n",
    "```{article-info}\n",
    ":author: Altay Sansal\n",
    ":date: \"{sub-ref}`today`\"\n",
    ":read-time: \"{sub-ref}`wordcount-minutes` min read\"\n",
    ":class-container: sd-p-0 sd-outline-muted sd-rounded-3 sd-font-weight-light\n",
    "```\n",
    "\n",
    "In this tutorial, we will use a public SEG-Y file located in Amazon Web Services' (AWS)\n",
    "Simple Storage Service (S3), also known as a cloud object store.\n",
    "\n",
    "This dataset, the Parihaka 3D full angle stack (4.7 GB per volume including full angle\n",
    "and near, mid, and far stacks), is provided by New Zealand Petroleum and Minerals (NZPM).\n",
    "Available information and data acquisition details are accessible via the\n",
    "[SEG Wiki][seg wiki], the [New Zealand GNS website], and the [NZPM data portal].\n",
    "\n",
    "[seg wiki]: https://wiki.seg.org/wiki/Parihaka-3D\n",
    "[nzpm data portal]: http://data.nzpam.govt.nz/GOLD/system/mainframe.asp\n",
    "[new zealand gns website]: http://www.gns.cri.nz\n",
    "\n",
    "```{important}\n",
    "For plotting, the notebook requires [Matplotlib](https://matplotlib.org/) as a dependency.\n",
    "Please install it before executing using `pip install matplotlib` or `conda install matplotlib`.\n",
    "```\n",
    "\n",
    "Let's start by importing some modules we will be using."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ff66db0926ba711",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from numpy.random import default_rng\n",
    "\n",
    "from segy import SegyFile\n",
    "from segy.config import SegyFileSettings\n",
    "from segy.schema import HeaderField\n",
    "from segy.standards import get_segy_standard"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc09a9cbef23c154",
   "metadata": {},
   "source": [
    "You can (but don't) download the SEG-Y directly clicking the [HTTP link] from the website.\n",
    "\n",
    "[http link]: http://s3.amazonaws.com/open.source.geoscience/open_data/newzealand/Taranaiki_Basin/PARIHAKA-3D/Parihaka_PSTM_full_angle.sgy\n",
    "\n",
    "This link is convenient as the `segy` library supports HTTP and we can directly use it\n",
    "without downloading as well. Hovewer, For demonstration purposes, we'll use the\n",
    "corresponding S3 link (or called bucket and prefix):\n",
    "\n",
    "`s3://open.source.geoscience/open_data/newzealand/Taranaiki_Basin/PARIHAKA-3D/Parihaka_PSTM_full_angle.sgy`\n",
    "\n",
    "It's important to note that the file isn't downloaded but rather read on demand from the\n",
    "S3 object store with the `segy` library.\n",
    "\n",
    "The `SegyFile` class uses information from the binary file header to construct a SEG-Y\n",
    "spec, allowing it to read the file. The SEG-Y Revision is inferred from the binary\n",
    "header by default, but can be manually set by providing a custom spec or adjusting settings.\n",
    "\n",
    "Since this is a public bucket and an object, we need to tell `S3` that we want anonymous\n",
    "access, which is done by configuring `storage_options` in settings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f83a983066a23eee",
   "metadata": {},
   "outputs": [],
   "source": [
    "path = \"s3://open.source.geoscience/open_data/newzealand/Taranaiki_Basin/PARIHAKA-3D/Parihaka_PSTM_full_angle.sgy\"\n",
    "\n",
    "# Alternatively via HTTP\n",
    "# path = \"http://s3.amazonaws.com/open.source.geoscience/open_data/newzealand/Taranaiki_Basin/PARIHAKA-3D/Parihaka_PSTM_full_angle.sgy\"\n",
    "\n",
    "settings = SegyFileSettings(storage_options={\"anon\": True})\n",
    "\n",
    "sgy = SegyFile(path, settings=settings)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba606aac90c33574",
   "metadata": {},
   "source": [
    "Let's investigate the JSON version of the SEG-Y spec for this file.\n",
    "\n",
    "Some things to note:\n",
    "1. The opening processed inferred Revision 1 from the binary header automatically.\n",
    "2. It generated the spec using **default** SEG-Y Revision 1 headers.\n",
    "3. We can check that some headers can be defined in the wrong byte locations.\n",
    "4. There are too many headers to deal with in the default schema.\n",
    "\n",
    "Note that we can build this JSON independently, load it into the spec\n",
    "and open any SEG-Y with a schema."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a24feb0fd760a8b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "sgy.spec"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64f7aa797249034f",
   "metadata": {},
   "source": [
    "Let's check the file size, number of traces, sample rate, etc. As expected, the file size\n",
    "matches what was in the description. We also observe there are ~ 1 million traces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61ee9a641245f155",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"file size: {sgy.file_size / 1024**3:0.2f} GiB\")\n",
    "print(f\"num traces: {sgy.num_traces:,}\")\n",
    "print(f\"sample rate: {sgy.sample_interval}\")\n",
    "print(f\"num samples: {sgy.samples_per_trace}\")\n",
    "print(f\"sample labels: {sgy.sample_labels // 1000}\")  # microsecond to millisecond"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3684405bdbfd2d7",
   "metadata": {},
   "source": [
    "Using the `SegyFile` we can read SEG-Y components.\n",
    "\n",
    "Here we read:\n",
    "- Textual file header\n",
    "- Binary file header\n",
    "- 1,000 traces (headers + data) from somewhere in the middle of the file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdf4c3cdd47e96e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "text_header = sgy.text_header\n",
    "binary_header = sgy.binary_header\n",
    "\n",
    "start = 500_000\n",
    "stop = start + 1_000\n",
    "\n",
    "traces = sgy.trace[start:stop]\n",
    "\n",
    "trace_headers = traces.header\n",
    "trace_data = traces.sample"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4cb6f58455ffa06",
   "metadata": {},
   "source": [
    "This should take around one second or less, based on internet connection.\n",
    "\n",
    "Let's print the textual header. There are a not many headers of interest. The\n",
    "available headers appear to be in the Revision 1 byte locations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb524f93a740fb90",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(text_header)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ba3bd7911a900ec",
   "metadata": {},
   "source": [
    "We can look at headers (by default it is a Pandas `DataFrame`) in a nicely formatted table.\n",
    "\n",
    "We can also do typical Pandas analytics (like plots, statistics, etc.) but it won't be shown here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f48a3e8ed773f45",
   "metadata": {},
   "outputs": [],
   "source": [
    "binary_header.to_dataframe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6cbdee064cef681f",
   "metadata": {},
   "outputs": [],
   "source": [
    "trace_headers.to_dataframe()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "368532bca3b3e876",
   "metadata": {},
   "source": [
    "Let's plot the traces.\n",
    "\n",
    "Note that they are all parsed from IBM floats to IEEE floats (decoded) in the background."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "237c21b1fbcc99d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12, 8))\n",
    "plot_kw = {\"aspect\": \"auto\", \"cmap\": \"gray_r\", \"interpolation\": \"bilinear\"}\n",
    "plt.imshow(trace_data.T, vmin=-1000, vmax=1000, **plot_kw);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f050ae1104240b38",
   "metadata": {},
   "source": [
    "## With Custom Spec\n",
    "\n",
    "We will create a new custom spec based on SEG-Y revision 1 with different binary\n",
    "and trace header fields. This way, we will parse ONLY the parts we want, with\n",
    "correct byte locations.\n",
    "\n",
    "A user can define a completely custom SEG-Y spec from scratch as well, but for\n",
    "convenience, we are customizing the Revision 1 schema with the parts that we\n",
    "want to change.\n",
    "\n",
    "Note that doing this will also modify the `segyStandard` field to `None` to\n",
    "ensure we don't assume the file schema is standard after doing this.\n",
    "\n",
    "From the binary file header, we will read:\n",
    "- Number of samples\n",
    "- Sample rate\n",
    "\n",
    "From the trace headers, we will read:\n",
    "- Inline\n",
    "- Crossline\n",
    "- CDP-X\n",
    "- CDP-Y\n",
    "- Coordinate scalar\n",
    "\n",
    "Based on the text header lines:\n",
    "\n",
    "```\n",
    "C 2 HEADER BYTE LOCATIONS AND TYPES:\n",
    "C 3     3D INLINE : 189-192 (4-BYTE INT)    3D CROSSLINE: 193-196 (4-BYTE INT)\n",
    "C 4     ENSEMBLE X: 181-184 (4-BYTE INT)    ENSEMBLE Y  : 185-188 (4-BYTE INT)\n",
    "```\n",
    "\n",
    "As we know by the SEG-Y Rev1 definition, the coordinate scalars are at byte 71."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec22a69697ebffd9",
   "metadata": {},
   "outputs": [],
   "source": [
    "rev1 = get_segy_standard(1.0)\n",
    "custom_spec = rev1.customize(\n",
    "    binary_header_fields=[\n",
    "        HeaderField(name=\"sample_interval\", byte=17, format=\"int16\"),\n",
    "        HeaderField(name=\"samples_per_trace\", byte=21, format=\"int16\"),\n",
    "        HeaderField(name=\"data_sample_format\", byte=25, format=\"int16\"),\n",
    "        HeaderField(name=\"num_extended_text_headers\", byte=305, format=\"int16\"),\n",
    "    ],\n",
    "    trace_header_fields=[\n",
    "        HeaderField(name=\"inline\", byte=189, format=\"int32\"),\n",
    "        HeaderField(name=\"crossline\", byte=193, format=\"int32\"),\n",
    "        HeaderField(name=\"cdp_x\", byte=181, format=\"int32\"),\n",
    "        HeaderField(name=\"cdp_y\", byte=185, format=\"int32\"),\n",
    "        HeaderField(name=\"coordinate_scalar\", byte=71, format=\"int16\"),\n",
    "    ],\n",
    ")\n",
    "\n",
    "sgy = SegyFile(path, spec=custom_spec, settings=settings)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de33a63c81b2072c",
   "metadata": {},
   "source": [
    "Now let's look at the spec again. It is a lot more compact."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca974e7209c6249",
   "metadata": {},
   "outputs": [],
   "source": [
    "sgy.spec"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "109579c0124ba864",
   "metadata": {},
   "source": [
    "As mentioned earlier, the JSON can be laded into the spec from a file too.\n",
    "\n",
    "```python\n",
    "from segy.schema import SegySpec\n",
    "import os\n",
    "\n",
    "json_path = \"...\"\n",
    "\n",
    "with open(json_path, mode=\"r\") as fp:\n",
    "    data = fp.read()\n",
    "    spec = SegySpec.model_validate_json(data)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e57eedd6f49eb172",
   "metadata": {},
   "source": [
    "Let's do something a little more interesting now. Let's try to plot a time slice\n",
    "by randomly sampling the file.\n",
    "\n",
    "We will read 5,000 random traces. This should take about 15-20 seconds, based on\n",
    "your internet connection speed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5d1855359930d69",
   "metadata": {},
   "outputs": [],
   "source": [
    "text_header = sgy.text_header\n",
    "binary_header = sgy.binary_header\n",
    "\n",
    "rng = default_rng()\n",
    "indices = rng.choice(sgy.num_traces, size=5_000, replace=False)\n",
    "\n",
    "traces = sgy.trace[indices]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1136a5e2f2143b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "binary_header.to_dataframe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f29d0c5c0595268",
   "metadata": {},
   "outputs": [],
   "source": [
    "trace_headers = traces.header.to_dataframe()\n",
    "\n",
    "trace_headers[\"cdp_x\"] /= trace_headers[\"coordinate_scalar\"].abs()\n",
    "trace_headers[\"cdp_y\"] /= trace_headers[\"coordinate_scalar\"].abs()\n",
    "\n",
    "trace_headers"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "686a354d6ede74c4",
   "metadata": {},
   "source": [
    "Now we can plot the time slice on actual coordinates from the headers and\n",
    "even see a hint of the outline of the data! Since we significantly\n",
    "down-sampled the data, the time slice is aliased and not very useful,\n",
    "but this shows us the concept of making maps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abc830a2a5f38c00",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 8))\n",
    "\n",
    "z_slice_index = 500\n",
    "\n",
    "x, y, z = (\n",
    "    trace_headers[\"cdp_x\"],\n",
    "    trace_headers[\"cdp_y\"],\n",
    "    traces.sample[:, z_slice_index],\n",
    ")\n",
    "\n",
    "scatter_kw = {\"ec\": [0.0, 0.0, 0.0, 0.5], \"linewidth\": 0.5}\n",
    "color_kw = {\"cmap\": \"gray_r\", \"vmin\": -1000, \"vmax\": 1000}\n",
    "\n",
    "plt.tripcolor(x, y, z, shading=\"gouraud\", **color_kw)\n",
    "plt.scatter(x, y, s=4, c=z, label=\"trace_locations\", **scatter_kw, **color_kw)\n",
    "plt.title(\"Water Bottom Map\")\n",
    "plt.colorbar(label=\"Amplitude\")\n",
    "plt.xlabel(\"X-Coordinate\")\n",
    "plt.ylabel(\"Y-Coordinate\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb01e8d1-efc1-474f-b73a-9b3a8867bbce",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}