ASAP_notebook_2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Anomaly Screening Applying Pseudo-wells (ASAP)\n",
    "\n",
    "Stavanger Machine Learning Hackathon 18-19 September, 2018.\n",
    "\n",
    "Project page: <https://events.agilescientific.com/project/asap>\n",
    "\n",
    "## Seismic data loading and processing\n",
    "\n",
    "This notebook shows how to:\n",
    "\n",
    "* (...)\n",
    "\n",
    "## To-do\n",
    "\n",
    "* (...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import segyio\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.interpolate import griddata\n",
    "%config InlineBackend.figure_format = 'retina'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import segypy\n",
    "\n",
    "# Set verbose level\n",
    "segypy.verbose=1;\n",
    "\n",
    "filename='./data/3d_nearstack.sgy';\n",
    "\n",
    "# Get only SegyHeader\n",
    "SH = segypy.getSegyHeader(filename);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inl=np.arange(1300,1502,2)\n",
    "crl=np.arange(1500,2002,2)\n",
    "hrz=np.recfromtxt('./data/Top_Heimdal_subset.txt', names=['il','xl','z'])\n",
    "\n",
    "filename='./data/3d_nearstack.sgy'\n",
    "#f=open(filename, 'rb')\n",
    "seis,header,trace_headers=segypy.readSegy(filename)\n",
    "near=seis.reshape(header['ns'],inl.size,crl.size)\n",
    "\n",
    "lagtime=trace_headers['LagTimeA'][0]*-1\n",
    "twt=np.arange(lagtime,header['dt']/1e3*header['ns']+lagtime,header['dt']/1e3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.interpolate import interp2d\n",
    "horizon = np.zeros((len(inl), len(crl)))\n",
    "print(horizon.shape)\n",
    "for i, idx in enumerate(inl):\n",
    "    for j, xdx in enumerate(crl):\n",
    "        time = hrz['z'][np.where((hrz['il']== idx) & (hrz['xl'] == xdx))]\n",
    "        if len(time) == 1:\n",
    "            horizon[i, j] = time \n",
    "            \n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(12, 12))\n",
    "ax.imshow(horizon, vmin=2030, vmax=2200)\n",
    "print(horizon)\n",
    "print(horizon[:, 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "points = []\n",
    "wanted = []\n",
    "for i in range(horizon.shape[0]):\n",
    "    for j in range(horizon.shape[1]):\n",
    "        if horizon[i, j] != 0.:\n",
    "            points.append([i, j, horizon[i, j]])\n",
    "        else:\n",
    "            wanted.append([i, j])\n",
    "points = np.array(points)\n",
    "print(points.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zs2 = scipy.interpolate.griddata(points[:, 0:2], points[:, 2], wanted, method=\"cubic\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for p, val in zip(wanted, zs2):\n",
    "    horizon[p[0], p[1]] = val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'horizon' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-22-c26c003fe637>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[0mfig\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0max\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msubplots\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfigsize\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m12\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m12\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0max\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mimshow\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mhorizon\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mvmin\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m2030\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mvmax\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m2200\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      3\u001b[0m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mhorizon\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mhorizon\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'horizon' is not defined"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABZUAAAVdCAYAAACPdqdzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzs3W/o7nddx/HXO//UXNvUcgV6w1zOHdASz/IPK1IHa3gjzCy64V/qhjpQRG8EC3SCBhI6ZgpRqCxDoUAEMxjpMnOJdUwqOdt0cWwxy9RwNudK9unGdR04/txv53V+v6vjER8P+PE51/fzvT7X59x98uXznbVWAAAAAACg8UPf6w0AAAAAAPD9Q1QGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAajuJyjPzwpl5x8x8Ymbunpk1M+874FqPm5l3z8xdM3PfzJyYmetn5lG72CsAAAAAAAf30B2t8ztJfjbJfyf5tySXHWSRmbkkyS1JLk7yoSS3Jnl6ktckuXpmrlhrfXUnOwYAAAAA4Izt6viL1ya5NMmFSV55iHXelU1QfvVa6/lrrd9eaz03yduTPCnJmw+9UwAAAAAADmzWWrtdcObZSW5O8idrrRedwfeekOSOJCeSXLLWuv+UuQuSfCnJJLl4rXXPLvcMAAAAAEDnXHpR33O3402nBuUkWWt9I8knkzwiyTPP9sYAAAAAANjY1ZnKu/Ck7Xj7PvOfT3JVNsdsfPTBFpqZY/tMPTmbc59PHGB/AAAAAADnkscnuXut9VNn80fPpah80Xb8+j7zJ68/8hC/8ZDzzjvv0UeOHHn0IdYAAAAAAPieO378eO69996z/rvnUlQ+ndmOpz0Eeq119AEXmDl25MiRpx07tt+DzAAAAAAA3x+OHj2az3zmMyfO9u+eS2cqn3wS+aJ95i/ccx8AAAAAAGfZuRSVb9uOl+4z/8TtuN+ZywAAAAAA/D87l6Lyzdvxqpn5jn3NzAVJrkhyb5JPne2NAQAAAACwcdaj8sw8bGYum5lLTr2+1rojyU3ZvLHwmj1fuy7J+UluXGvdc1Y2CgAAAADAd9nJi/pm5vlJnr/9+JPb8Vkz897tv7+y1nr99t+PTXI8yRezCcinelWSW5LcMDNXbu97RpLnZHPsxbW72C8AAAAAAAezk6ic5KlJXrrn2hO2f8kmIL8+p7HWumNmLk/ypiRXJ3leki8luSHJdWutr+1ovwAAAAAAHMBOovJa641J3ljeeyLJPMj8nUlevot9AQAAAACwW+fSi/oAAAAAADjHicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABAbWdReWYeNzPvnpm7Zua+mTkxM9fPzKPOcJ2fn5kPbb//rZn515n5yMxcvau9AgAAAABwMDuJyjNzSZJjSV6e5NNJ3p7kX5K8JsnfzsyPleu8Msknkly5Hd+e5ONJfjHJX8zMtbvYLwAAAAAAB/PQHa3zriQXJ3n1WusdJy/OzNuSvDbJm5O84sEWmJmHJfndJN9KcnStddspc29J8g9Jrp2Z31tr3bejfQMAAAAAcAYO/aTyzDwhyVVJTiR5557pNyS5J8mLZ+b80yz16CQXJbn91KCcJGut40luT3Jekh897J4BAAAAADiYXRx/8dzteNNa6/5TJ9Za30jyySSPSPLM06zz5ST/meTSmXniqRMzc2mSJyb57FrrqzvYMwAAAAAAB7CL4y+etB1v32f+89k8yXxpko/ut8haa83MNUnel+TYzHwwyV1JHpvkV5J8LslvNBuamWP7TF3WfB8AAAAAgAe2i6h80Xb8+j7zJ68/8nQLrbX+dGbuSvL+JC85Zeo/krwnm5f/AQAAAADwPbKL4y9OZ7bjOu2NMy9K8pdJPpHkSDbHZhzJ5gnn30/ygeYH11pHH+gvya0H+Q8AAAAAALCxi6h88knki/aZv3DPfQ9oe27yu7M55uLFa61b11r3rrVuTfLiJMeS/NrMPPvwWwYAAAAA4CB2EZVv246X7jN/8qV7+525fNJVSR6W5OMP8MK/+5P89fbj0YNsEgAAAACAw9tFVL55O141M9+x3sxckOSKJPcm+dRp1vnh7fiYfeZPXv+fg2wSAAAAAIDDO3RUXmvdkeSmJI9Pcs2e6euSnJ/kxrXWPScvzsxlM3PZnns/sR1fODM/c+rEzDw1yQuzOZf5Y4fdMwAAAAAAB/PQHa3zqiS3JLlhZq5McjzJM5I8J5tjL67dc//x7XjyJX5Za316Zt6T5OVJ/m5mPpjki9nE6ucneXiS69dan9vRngEAAAAAOEM7icprrTtm5vIkb0pydZLnJflSkhuSXLfW+lq51G9mc3byy5L8UpILktyd5G+S/OFa6wO72C8AAAAAAAezqyeVs9a6M5unjJt7Z5/rK8l7t38AAAAAAJxjdvGiPgAAAAAAfkCIygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBtZ1F5Zh43M++embtm5r6ZOTEz18/Mow6w1lNm5saZuXO71pdn5uMz85Jd7RcAAAAAgDP30F0sMjOXJLklycVJPpTk1iRPT/KaJFfPzBVrra+Wa70syR8l+WaSDyc5keSRSZ6c5HlJbtzFngEAAAAAOHM7icpJ3pVNUH71WusdJy/OzNuSvDbJm5O84nSLzMwzswnK/5zk6rXWv++Zf9iO9gsAAAAAwAEc+viLmXlCkquyeaL4nXum35DkniQvnpnzi+XemuQhSV60NygnyVrrfw+3WwAAAAAADmMXTyo/dzvetNa6/9SJtdY3ZuaT2UTnZyb56H6LzMzjkvxCkr9P8rmZeU6So0lWks8muXnv+gAAAAAAnF27iMpP2o637zP/+Wyi8qV5kKic5OdOuf9jSZ69Z/6fZuYFa60vnG5DM3Nsn6nLTvddAAAAAAD2d+jjL5JctB2/vs/8yeuPPM06F2/HX09yJMkLtmv/dJI/TvKUJH8+Mw8/+FYBAAAAADiMXb2o78HMdlynue8hp4y/tdb68Pbz3TPz0mxC8+VJfjXJ+x9sobXW0QfcyOYJ5qc1mwYAAAAA4Lvt4knlk08iX7TP/IV77tvPf23H+5J85NSJtdZK8qHtx6ef6QYBAAAAANiNXUTl27bjpfvMP3E77nfm8t51vrHPC/lORufzzmBvAAAAAADs0C6i8s3b8aqZ+Y71ZuaCJFckuTfJp06zzj8m+UqSH5+Zn3iA+SdvxxMH3yoAAAAAAIdx6Ki81rojyU1JHp/kmj3T1yU5P8mNa617Tl6cmctm5rI963w7yR9sP7711EA9M09J8rIk307yZ4fdMwAAAAAAB7OrF/W9KsktSW6YmSuTHE/yjCTPyebYi2v33H98O86e629JcmWSlyR5ysz8VZLHZPNyvh9J8rq11hd2tGcAAAAAAM7QLo6/OPm08uVJ3ptNTH5dkkuS3JDkWWutr5brfDObqHxdkkdk8+TzL2cTrJ+31nrbLvYLAAAAAMDB7OpJ5ay17kzy8vLevU8onzr3zSRv3P4BAAAAAHAO2cmTygAAAAAA/GAQlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAID/Y+9+Qj676juOf75Fusg0jrrISqgYksw6ilqDqUlgCIIoFinUiiYrSUpE3AgRJKXipn+mEYUqTFNrW2k2cVGQgMSNwYWD2SUxRGKFtLZElDC1Vs1x8fwenD76ZD7z/G6blrxe8OPAPfd37lm/OdwLAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAFr1csAAAgAElEQVQAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQE5UBAAAAAKiJygAAAAAA1ERlAAAAAABqojIAAAAAADVRGQAAAACAmqgMAAAAAEBNVAYAAAAAoCYqAwAAAABQ2ywqz8xrZ+b8zDw7Mz+ZmWdm5tzMvHqPNW+emZ/PzJqZP9lqrwAAAAAAnMwrtlhkZq5N8miSa5J8OckTSd6U5MNJbp+Zm9Zaz13hmlcn+Zsk/5Hkt7bYJwAAAAAA+9nqpPJncxCU71lrvXut9bG11q1J/iLJDUk+eYI1/zLJ6SSf2miPAAAAAADsae+oPDOvT3I2yTNJPnNk+hNJLiZ5/8ycuoI135XkjiT3JHl23z0CAAAAALCNLU4q37obH15rvXDpxFrr+SRfT3JVkrc0i83MNUk+n+ShtdYXN9gfAAAAAAAb2eKdyjfsxm8fM/9UDk4yX5/kq8V6n8tB7P7QSTc0MxeOmTpz0jUBAAAAANgmKp/ejT86Zv7w+qsut9DM3JnkXUl+f631/Q32BgAAAADAhraIypczu3G96E0zr0tyLsmDa61/3OeBa603HPOMC0lu3GdtAAAAAICXsy3eqXx4Evn0MfOvPHLfcc4n+XGSuzbYEwAAAAAA/wO2iMpP7sbrj5m/bjce987lQzcmuSbJv8/MOvwl+evd/L27aw/tt10AAAAAAE5qi9dfPLIbz87Mb6y1XjicmJmrk9yUgxPI37jMOl9IctWvuX5dkpuTPJbkQpJv7b1jAAAAAABOZO+ovNZ6emYeTnI2yd1JPn3J9H1JTiX5q7XWxcOLM3Nm998nLlnnnl+3/sx8MAdR+Z/WWh/fd78AAAAAAJzcVh/quyvJo0nun5nbkjye5M1JbsnBay/uPXL/47txAgAAAADA/xtbvFM5a62nk7wxyQM5iMkfTXJtkvuT/M5a67ktngMAAAAAwEtrq5PKWWt9L8kd5b31CeW11gM5iNUAAAAAALzENjmpDAAAAADAy4OoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANREZQAAAAAAaqIyAAAAAAA1URkAAAAAgJqoDAAAAABATVQGAAAAAKAmKgMAAAAAUBOVAQAAAACoicoAAAAAANQ2i8oz89qZOT8zz87MT2bmmZk5NzOvLv9/ambeNzN/PzNPzMzFmXl+Zr45Mx+dmd/caq8AAAAAAJzMK7ZYZGauTfJokmuSfDnJE0nelOTDSW6fmZvWWs9dZpm3Jflikh8keSTJQ0lek+SdSf40yXtm5ra11n9usWcAAAAAAK7cJlE5yWdzEJTvWWt9+vDizPx5ko8k+WSSD11mjX9N8odJHlxr/dcla1yd5GtJ3prk7iR/ttGeAQAAAAC4Qnu//mJmXp/kbJJnknzmyPQnklxM8v6ZOfVi66y1Hltr/d2lQXl3/fn8MiS/fd/9AgAAAABwclu8U/nW3fjwWuuFSyd2QfjrSa5K8pY9nvHT3fizPdYAAAAAAGBPW7z+4obd+O1j5p/KwUnm65N89YTPuHM3fqW5eWYuHDN15oTPBwAAAAAg25xUPr0bf3TM/OH1V51k8Zn5oyS3J3ksyfmTrAEAAAAAwDa2+lDfi5nduK74jzPvSXIuBx/x+7211k8v85eDB631hmPWu5DkxivdBwAAAAAAB7Y4qXx4Evn0MfOvPHJfZWbeneRLSf4tydvXWt852fYAAAAAANjKFlH5yd14/THz1+3G4965/Ctm5r1JHkzy/SS/u9Z68jJ/AQAAAADgf8EWUfmR3Xh2Zv7bejNzdZKbkvw4yTeaxWbmD5L8Q5JncxCUn9pgjwAAAAAAbGDvqLzWejrJw0lel+TuI9P3JTmV5AtrrYuHF2fmzMycObrWzHwgyd8m+eckN3vlBQAAAADA/y1bfajvriSPJrl/Zm5L8niSNye5JQevvbj3yP2P78bDj/hlZm5Jcj4HofuRJHfMzJG/5YdrrXMb7RkAAAAAgCu0SVReaz09M29M8sdJbk/yjiT/kuT+JPettX5QLPPb+eXJ6TuPuee7SURlAAAAAICXyFYnlbPW+l6SO8p7f+UI8lrrgSQPbLUfAAAAAAC2t8WH+gAAAAAAeJkQlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAQE1UBgAAAACgJioDAAAAAFATlQEAAAAAqInKAAAAAADURGUAAAAAAGqiMgAAAAAANVEZAAAAAICaqAwAAAAAwC/au/9YS8ryDuDfh4IGcOWnKAkmVMqypFKIpAhSkYUESRrDRqA1VgSNaVolNCppmkoETGmTtgqtPxJrs02xtqQ0ARKVQiqoFEJUkBLNggpdlYKi0MqPAi3u2z9mbj1euLuz9845957l80luJufMnGfe2c19Mud7Z94ZTKgMAAAAAMBgQmUAAAAAAAYTKgMAAAAAMJhQGQAAAACAwYTKAAAAAAAMJlQGAAAAAGAwoTIAAAAAAIMJlQEAAAAAGEyoDAAAAADAYEJlAAAAAAAGEyoDAAAAADCYUBkAAAAAgMGEygAAAAAADCZUBgAAAABgMKEyAAAAAACDCZUBAAAAABhMqAwAAAAAwGBCZQAAAAAABhMqAwAAAAAwmFAZAAAAAIDBhMoAAAAAAAwmVAYAAAAAYDChMgAAAAAAgwmVAQAAAAAYTKgMAAAAAMBgQmUAAAAAAAYTKgMAAAAAMJhQGQAAAACAwYTKAAAAAAAMJlQGAAAAAGAwoTIAAAAAAIONFipX1SFVtbmqHqyqZ6pqa1VdUVX77WSd/fvPbe3rPNjXPWSssQIAAAAAsDy7j1Gkqg5LcluSg5Jcl+SeJMcl+b0kp1fVia21RwbUOaCvsz7JTUmuSrIhyTuS/HpVndBau3+MMQMAAAAAsPPGulL5E+kC5Qtaa5taa3/QWjslyeVJjkhy2cA6f5wuUL68tXZqX2dTunD6oH4/AAAAAACskhWHylX1qiSnJdma5OOLVl+c5Mkk51TV3juos3eSc/rtL160+mN9/Tf2+wMAAAAAYBWMcaXyKf3yxtbatskVrbXHk9yaZK8kx++gzglJ9kxya/+5yTrbktzYv9y44hEDAAAAALAsY8ypfES//NYS67+d7krm9Um+sMI66etsV1XdscSqo7ds2ZJjjz12RyUAAAAAANa0LVu2JMmhs97vGKHyPv3yJ0usX3h/3xnV2Z7dnnrqqZ/eeeed/7aCGsALw4Z+ec+qjgKYB/oFMJR+AQylXwBDHZ3kJbPe6Rih8o5Uv2yzqtNae95LkReuYF5qPcAC/QIYSr8AhtIvgKH0C2Co7czYMFVjzKm8cAXxPkusf+mi7aZdBwAAAACAKRkjVL63Xy411/Hh/XKpuZLHrgMAAAAAwJSMESrf3C9Pq6qfq1dV65KcmOSpJLfvoM7t/XYn9p+brLNbuof9Te4PAAAAAIAZW3Go3Fq7L8mN6Z4y+J5Fqy9NsneSK1trTy68WVUbqmrD5IattSeSfLrf/pJFdc7v69/QWrt/pWMGAAAAAGB5xnpQ37uT3JbkL6vq1CRbkrw2ycZ001V8YNH2W/plLXr/D5OcnOR9VXVMkq8kOTLJGUkeznNDawAAAAAAZqhaa+MUqnplkg8lOT3JAUkeSnJtkktba48u2rYlSWttcaicqto/ycVJNiU5OMkjSa5P8sHW2gOjDBYAAAAAgGUZLVQGAAAAAGDXN8aD+gAAAAAAeIEQKgMAAAAAMJhQGQAAAACAwYTKAAAAAAAMJlQGAAAAAGAwoTIAAAAAAIMJlQEAAAAAGGzuQ+WqOqSqNlfVg1X1TFVtraorqmq/nayzf/+5rX2dB/u6h0xr7MBsrbRfVNXeVfVbVfX3VXVPVT1ZVY9X1deq6v1V9aJpHwMwG2OdXyyqeVJV/bSqWlX90ZjjBVbHmL2iqo6qqiur6vt9rYer6ktV9fZpjB2YrRGzi1+rquv6zz9dVd+rqs9X1enTGjswO1V1VlV9tKpuqarH+u8Of7fMWqN/p/m5+q21Meqsiqo6LMltSQ5Kcl2Se5Icl2RjknuTnNhae2RAnQP6OuuT3JTkq0k2JDkjycNJTmit3T+NYwBmY4x+0Z+oXZ/k0SQ3J/lOkv2TvCnJK/r6p7bWnp7SYQAzMNb5xaKa65LcneTAJC9Jcllr7aIxxw3M1pi9oqrOS/LXSf47yWeTbE2yb5JXJ3mwtfaWkYcPzNCI2cXvJvlEkieTXJPkgSSHJHlzkr2SXNRau2waxwDMRlXdleToJE+k+x3fkOQzrbW37WSd0b/TPGcfcx4q35DktCQXtNY+OvH+R5K8N8knW2u/M6DOJ5P8dpLLW2vvm3j/giR/keSG1pq/+sEcG6NfVNUxSX45ydWttf+ZeH9dki8meU2SC1trHx7/CIBZGev8YlHNzUk2JfnzJJdFqAxzb8TvIscn+dck30hyemvtB4vW79Fa+99RBw/M1EjfRfZI8qMkL05yTGvt3ol1Ryb5epJtSfZrrT0z/lEAs1BVG9OFyd9J8oZ0F7QtJ1Qe/TvNc/Yxr6FyVb0qyX3p/op/WGtt28S6dUkeSlJJDmqtPbmdOnuna8zbkhzcWnt8Yt1u/T4O7ffhamWYQ2P1ix3s461JPpPks621N6140MCqmEa/qKozklyb5Jwkuyf5mwiVYa6N2Suq6stJXp/kqNbaN6Y2aGBVjJhdvDzJD5Lc3Vo7+nnW353kqCQHrvTqQ2BtqKqTs4xQeRYZSDLfcyqf0i9vnPzHSZI+GL413e0fx++gzglJ9kxy62Sg3NfZluTG/uXGFY8YWC1j9YvtWbiC6NkV1ABW36j9oqoOSvKpJNe21pY1FxqwJo3SK/rnt7w+ydeSfLOqNlbVhf2zGk7tL3IB5ttY5xYPp7sgbn1VHT65oqrWJzk8yV0CZSCzyUDmOlQ+ol9+a4n13+6X62dUB1i7ZvF7/s5++c8rqAGsvrH7xV+lO99a0a1lwJozVq/41Yntb+p//izdVDn/kuSuqvqlFYwTWH2j9IvW3Wb+nnTnFXdU1d9W1Z9U1ZVJ7kjyzSRnjzBeYP7NJOvcfSUfXmX79MufLLF+4f19Z1QHWLum+nteVecnOT3JXUk2L6cGsGaM1i+q6p3pHvr7m621H44wNmDtGKtXHNQvfyPJj9M9bOsLSV6W5OJ00+Z8rqqOmnyeAzBXRju3aK1dXVUPJvmHJG+fWPXDdNNrmbITSGaUdc7zlco7Uv1ypZNGj1UHWLuW/XteVW9OckW6+c3O9CAd2OUN6hdVdWi63nB1a+0fpzwmYO0Zem7xCxPLd7XWrmmtPdZauy/JuemmxVif5MzpDBNYAwZ/F6mqt6W7i+GWJEemu339yHR/jPpYkqumNEZg1zJK1jnPofJCqr7PEutfumi7adcB1q6p/J5X1aZ0J24PJznZwzxhlzBWv9ic5Kkk7x5jUMCaM1av+M9++UySz0+u6G91v65/edzODhBYM0bpF/28yZvTTXNxTmvtntbaU621e9Ld1XBHkrP7B3sBL2wzyTrnOVS+t18uNf/HwsT1S80fMnYdYO0a/fe8qs5OcnW6W83e0Fq7dwcfAebDWP3iNelua/9RVbWFn3S3pibJB/r3rl3ZcIFVMvZ3kccXP0intxA677kTYwPWlrH6xWlJ9kjyped58Na2JF/uXx67nEECu5SZZJ3zPKfyzf3ytKrabbKpVtW6JCemu0Lo9h3Uub3f7sSqWtc/BXGhzm7pGvfk/oD5M1a/WPjMW5NcmeQ/kmx0hTLsUsbqF1emuyV1scOTnJRuDvY7knx9xSMGVsNYveLudHMpH1hVL3+e+ddf3S+3rnzIwCoZq1+8uF++bIn1C++bfx0YNQNZytxeqdzPM3ZjkkPTPQF10qVJ9k5yZWvtyYU3q2pDVW1YVOeJJJ/ut79kUZ3z+/o3CI1gfo3VL/r3z03XM76X5CS9AXYtI55fXNBae9fin/zsSuXP9e99fGoHA0zNiL3i2SSf7F/+aX9Ry8L2RyU5L8mzSf5p5EMAZmTE7yK39MuzqupXJldU1TFJzko3P+pN440eWMuqao++Xxw2+f5y+s6y9t9N1TWf+n+029LdXnpdki1JXptkY7pLuF/XWntkYvuWJK21WlTngL7O+nQN+CvpJrs/I91cqa/r/0OAOTVGv6iqjekejOTCNqgAAAHESURBVLFbuvnMvv88u/qv1toVUzoMYAbGOr9YovZ56YLly1prF40+eGBmRvwusle6h2wdn+7uhS+mu+LwzHTTXry/tfaRKR8OMEUj9ovNSd6R7mrka5J8N11otCnJi5Jc0Vp775QPB5ii/tlNm/qXr0jyxiT352d/WPpxa+3CfttDk/x7ku+21g5dVGen+s6yxjrPoXKSVNUrk3woyelJDkjyUJJrk1zaWnt00bZLfumrqv2TXJzuP+7gJI8kuT7JB1trD0zzGIDZWGm/mAiDtuc5zRyYP2OdXzxP3fMiVIZdxojfRfZK8vtJ3pLkF5M8neSrST7cWrt+mscAzMYY/aKqKsm56e5iODrJuiSPpfuD1Kdaa1dN9yiAaauqS9Llk0v5/8xhe6Fyv35w31nWWOc9VAYAAAAAYHbmdk5lAAAAAABmT6gMAAAAAMBgQmUAAAAAAAYTKgMAAAAAMJhQGQAAAACAwYTKAAAAAAAMJlQGAAAAAGAwoTIAAAAAAIMJlQEAAAAAGEyoDAAAAADAYEJlAAAAAAAGEyoDAAAAADCYUBkAAAAAgMGEygAAAAAADCZUBgAAAABgMKEyAAAAAACD/R8ZDlgMyMNEJgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x864 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 686,
       "width": 714
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(12, 12))\n",
    "ax.imshow(horizon, vmin=2030, vmax=2200)\n",
    "print(horizon)\n",
    "print(horizon[:, 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 38\n",
    "M = 0\n",
    "hrz_idx = [np.abs(twt-val).argmin() for val in horizon[N, :]]\n",
    "print(hrz_idx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(12, 12))\n",
    "top_add, below_add = 12, 50\n",
    "amplitude = near[hrz_idx[M]-top_add:hrz_idx[M]+below_add, N, M]\n",
    "ax.plot(amplitude, range(len(amplitude)))\n",
    "ax.axhline(top_add)\n",
    "ax.set_ylim(below_add+top_add, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h_bin = np.zeros((250, near[:, N].shape[1]))\n",
    "print(h_bin.shape)\n",
    "for i, val in enumerate(hrz_idx):\n",
    "    h_bin[val, i] = 1\n",
    "print(near.shape)\n",
    "fig, ax = plt.subplots(1, 1, figsize=(12, 12))\n",
    "clip = abs(np.percentile(seis, 0.8))\n",
    "ax.imshow(near[:, N], cmap=\"Greys\", vmin=-clip, vmax=clip)\n",
    "#ax.imshow(h_bin, alpha=0.2)\n",
    "ax.plot(range(len(hrz_idx)), hrz_idx, linewidth=5, color=\"black\")\n",
    "ax.axvline(276//2, color=\"red\", linewidth=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_add, below_add = 12, 52\n",
    "traces = np.zeros((horizon.shape[0], horizon.shape[1], top_add+below_add))\n",
    "for i in range(horizon.shape[0]):\n",
    "    hrz_idx = [np.abs(twt-val).argmin() for val in horizon[i, :]]\n",
    "    for j in range(horizon.shape[1]):\n",
    "        amplitude = near[hrz_idx[j]-top_add:hrz_idx[j]+below_add, i, j]\n",
    "        traces[i, j, :] = amplitude"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axarr = plt.subplots(2, 1, figsize=(12, 12))\n",
    "axarr[0].imshow(traces[0].T)\n",
    "axarr[1].imshow(traces[:, 0].T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first read in well data from Per Avseth's [open dataset distributed with the first edition of QSI (Quantitative Seismic Interpretation)](https://srb.stanford.edu/quantitative-seismic-interpretation).\n",
    "\n",
    "Well 2 corresponds to the Glitne well 15/5-5 drilled in 1995. Well 5 is 15/5-6 (1997) and well 3 is 15/5-3. For all the informations related to these wells check [NPD's Factpages](http://factpages.npd.no/factpages/)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}