{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I've noticed in some places that the model bathymetry in Boundary Pass doesn't seem to match with reality. \n",
    "\n",
    "See these two notebooks:\n",
    "* http://nbviewer.ipython.org/urls/bitbucket.org/salishsea/analysis/raw/tip/Nancy/strat/comparisons/ONC%20Moorings.ipynb\n",
    "* http://nbviewer.ipython.org/urls/bitbucket.org/salishsea/analysis/raw/tip/Nancy/strat/comparisons/Fall%202014%20IOS%20data.ipynb\n",
    "\n",
    "Is this a consistent error or is it just a discrepancy with the closest model grid point? I will investigate by looking at the bathymetry near the observations points as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "import pandas as pd\n",
    "sys.path.insert(0, './strat/comparisons/')\n",
    "import ONC\n",
    "import ACTDR\n",
    "import comparisons\n",
    "from salishsea_tools import tidetools\n",
    "import netCDF4 as nc\n",
    "import matplotlib.pyplot as plt\n",
    "import datetime\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "grid = nc.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc')\n",
    "\n",
    "bathy = grid.variables['Bathymetry'][:]\n",
    "lons = grid.variables['nav_lon'][:]\n",
    "lats = grid.variables['nav_lat'][:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ONC Mooring Point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "csvfilename = ('/ocean/nsoontie/MEOPAR/ONC/BoundaryPass//'\n",
    "               'BoundaryPass_BoundaryPassMooring_CTD_20140411T175025Z_20150408T235058Z-NaN_clean.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(48.766222, -123.039194, 234.0)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "./strat/comparisons/ONC.py:37: FutureWarning: convert_objects is deprecated.  Use the data-type specific converters pd.to_datetime, pd.to_timedelta and pd.to_numeric.\n",
      "  data = data.convert_objects(convert_numeric=True)\n"
     ]
    }
   ],
   "source": [
    "data, lat, lon, depth =  ONC.load_mooring_csv(csvfilename)\n",
    "print(lat,lon, depth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(343, 289, 205.0)\n"
     ]
    }
   ],
   "source": [
    "j,i = tidetools.find_closest_model_point(lon,lat,lons, lats, bathy)\n",
    "print(j,i, bathy[j,i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At the closest grid point, the model bathymetry is 30m shallower. What about neighbouring points?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[184.0 192.0 188.0 186.0 101.625 62.375]\n",
      " [197.0 196.0 208.0 208.0 178.0 116.625]\n",
      " [207.0 203.0 208.0 204.0 189.0 159.0]\n",
      " [202.0 205.0 210.0 205.0 180.0 201.0]\n",
      " [195.0 196.0 192.0 192.0 200.0 197.0]\n",
      " [195.0 191.0 186.0 139.25 164.625 192.0]]\n"
     ]
    }
   ],
   "source": [
    "stp = 3\n",
    "print(bathy[j-stp:j+stp,i-stp:i+stp])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f2bfc941410>]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/nsoontie/anaconda3/envs/py2/lib/python2.7/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison\n",
      "  if self._edgecolors == str('face'):\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAETCAYAAAAmkv2xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXvcHFWVrp+XJEC4KCIKBKIBCRoCSMCBIChBhIk4clEZ\nQVFBEc5IBNFRBDwnX9SZ8TIqgygGRQS5qYAZ5CI3ExRBhJCQQBIhSJQACYgGEy4hl3X+2NVJfZ3u\nruru6u6q6vXkV7907dq7alel89bqtddeW2aG4ziOU3w26nUHHMdxnGxwQXccxykJLuiO4zglwQXd\ncRynJLigO47jlAQXdMdxnJJQCEGX9FlJayVtXef4WZIekjRX0hWSNonKr5I0K9oekzQrxbV+Jenv\nkn6Z9X04juN0ktwIuqQJki6uUT4SOBT4c512o4BPAHub2R7AEOBYADM71szGmdk44JpoS+LrwIdb\nuQfHcZxekhtBB+rNcPoW8PkG7f4BrAI2kzQU2Ax4Il5BkoB/Ba6M9odI+oakP0h6QNLJ6zph9mtg\nReu34TiO0xvyJOjaoEA6ElhsZnPqNTKzvwHfBP4CPAksM7Pbqqq9DVhqZo9G+x+P6u0L7At8IrL0\nHcdxCsvQXndA0u+BTYAtgK1jfu7JwFnAYfHqNdq/Afg0MAp4Dvi5pA+Z2eWxascBV8T2DwP2kPT+\naP8VwC7Aonbvx3Ecp1f0XNDNbDyApIOAE8zsxGh/d2An4IHgMWFHYKakfc3s6dgp3gLcZWbPRu2u\nBd4KXB7tDwWOBvauuvQkM7u1XreyuDfHcZxukuhykTRR0gJJj0g6s8bxN0m6W9JLkj7bTNvqU8V3\nzOxBM9vWzHYys52AxYSBz6er2i0AxksaHvnK3wnMix1/JzDfzJ6Mld0MfDISeyTtKmmzen1xHMeB\nEKQhaXoUVfegpNOi8i9H43GzJd0eBXPE271O0opqjcyahoIuaQhwPjAR2A04TtKYqmrPAp8C/ruF\ntnGMxpbxumOSRki6AcDMHgAuBe4DKr72C2PtPkA0GBrjhwTRv1/SXOACQnQMkn4L/Aw4RNLjkg5t\n0CfHcfqLVcAZZjYWGA+cGuna183szWa2FzCN4DKO8y3ghk53Lsnlsi+w0MwWQYjrBo4E5lcqmNkz\nwDOS3t1s2zhmdgdwR72OmNnOsc9PAu+O7X+dEG5Yq92JNcoMOCfaqo+9rV4fHMfpb8xsCbAk+rxC\n0nxghJnFdW0L4K+VHUlHAX8Cnu90/5JcLjsAj8f2F0dlaWinreM4Tq6JIuPGAfdE+/8h6S/AR4Gv\nRmVbEMKuB7rRpyRBb2dw0AcWHccpJZFQXw2cbmYrAMzsHDN7HfBj4NtR1QHg22b2Al0Ym0tyuTwB\nxJ37IwmWdhpStZXkwu84TmrMrC1hbFZzqq8naRhh1vllZjatRpMrgBujz/sC75P0dWArYK2kF83s\ne833PJkkQb8PGB39tHiSMMB4XJ261Q85dVtfBi8wMDDAwMBAZueTpmR2rjSYVY8DbUjaPk2ebG09\ni2bvvVHf4+eqV6+zz3o6cHAHz18kBjI5y1dS1vti1X4USXcRMM/Mzo2VjzazR6LdI4FZAGb29lid\nycDyTok5JAi6ma2WNIkQ5jcEuMjM5ks6JTo+VdJ2wL2EyTlrJZ0O7BYNGGzQtlM34myI2eS6QhMX\npnbFKI2QN3uubop52nM1us9Gz9rJH8Nab3oAcDwwJzYJ8mzg45LeCKwBHgX+rc0utkTixCIzuwm4\nqapsauzzEga7Vhq2dXpPJwS4mfppXjKt0I6gSlMSX3LVdVrBhT8ftDqj0szupPbYY6LOmVnH/+F7\nPlPUWc+ECRMyP2dcQOqJUSsi0+mXQjPPolMWebN07mU1qo22Ti2G97oDHUK99l9Lsl73wUkvZFkJ\neRbWbuU83Sap3733uRcXu3Zg3We9d6BGjYFMBkV/kLLuJ2h/ELabuIXuAOms9F66aqopgiBm+bz6\ngbiYd5qyCp9b6M46OuXbzpK8CHlRf13kkUZCvqGVno2FflnKusfjFrpTUKqt9DwJObgAlo1uWuTV\nlFX48rTAhZMj8iTm0pTciXkW/cnTM+4FtX3kg+mU6A9LuRWNsr6onBbpd5Fx+oMiinUa3EJ3ck9e\nXzJupbdPr6z04Sm3ouEWuuM4fUdZhc8tdKcQ5NWSdSu9fdJY6X+yqXw9rDKZCe5DdxzH6TKPHb19\nR85bVuHzOHSnUOQt2qWCx6W3T7WvvJaY76ynMolDn5my7j54HLrj9B1ZpTJwOmeVxymr8LkP3SkU\nZRbNMt9bGvTega6IOZTXh+6C7jgZ0e8ukyLRatiipJGSpkt6SNKDkk6Lyr8s6QFJsyXdLmlkVL51\nVH+5pO90+r7ch+4UkryKp/vS2+dP65db2ICsfOhp19HckcE+9GhBn+3MbHa0ruhM4ChgsZktj+p8\nCnizmZ0kaTPCQtK7A7ub2afa6XsSbqE7hSSP7ok89smpzdCUWzVmtsTMZkefVwDzgREVMY/YAvhr\nVOcFM/sdsLJDtzKIso4NOE7XyFrI+31Vo511SkMrPQuGpVW+1fUPResljwPuifb/A/gw8AIwvqp6\nV9wQbqE7haXXFrHZ5J73wWmNoUPTbfWI3C1XA6dHljpmdo6ZvQ74MfDtLtzGBriF7jgt0Gkhdyu9\ns1b6sCG1y3+zJmyNkDQMuAa4zMym1ahyBXBjez1sDRd0p9B0W/jcIu8dO//6SSCbOT71rO93DIV3\nxPb/87nBxyUJuAiYZ2bnxspHm9kj0e6RwKyqU3dlcpJHuTiFpxuC3ish72crHULESxDyiEOUSZSL\nvTZl3ac3iHI5EPgNMIf1fvGzgY8DbwTWAI8C/2ZmT0dtFgFbAhsDfwcOM7MF7dxDPdxCdwpPp610\nt8p7xyAxz5IWlc/M7qT22ONNDdqMau1qzeOC7jh1yIOQ960vfcvo2R8F1PJSt0tJla+kt+X0G1kK\nXx6EvG/ZskvPvqTK52GLjhMjj2Kexz51hHpiflQHrjUk5VYwSvqecvqRdqz0vhHNPNItqzxOSZXP\nLXSnrynK5KAi9LFllqd4CWdtpW+ScisYLuhOqWhG+Eotkk5jWk3mknNc0J2+oyhWeTVF7HNqum2l\nu6A7TjFoJHylFkUnPSUdFHVBd/qColrl1ZThHuqSxkrPCrfQHac4VIQvCyHPcmJPPicJNcgRmxeW\nT8lW8Esq6AXssuOkI0shb3cR6HwK+TLgSoIMfKLHfYlYPmVwGGOnrPaSKl9Jb8txWqfTFnmrL4fs\n+rUauBu4C9gbuA/4B/CKjM6fEZ10wRQwJDENLuiOE6ORaDYrxPl01TwG3AC8imCVb01IGvhb4N0Z\nXaNNuuFLL6nyJd6WpInAuYQx3x+a2ddq1DkPeBdh6aUTzGxWVH4WcDywFpgLnGhmXVlbz3Gaodvi\nm/blkF2/lgO3AH8h/Fd9I+tTdB8AnB/9vVVG18s5BYxgSUPDQVFJQwj/0hOB3YDjJI2pqnM4sIuZ\njQZOBi6IykcRTIC9zWwPwiM8NuP+O05bSFOaEs0kCz4rAc7uXGuA3xP+W74SOBV4E4PXW9gc2IeQ\n5rtPaHFQVNJISdMlPSTpQUmnReXfkDRf0gOSrpX0yqh8U0lXSpojaZ6kL3TytpKiXPYFFprZIjNb\nBVxFWI0jzhHAJQBmdg+wlaRtCU65VcBmkoYCmwFPZNl5x2mHXrtE6rXJrl+PAz8A/gicCLyTsMZC\nLd5KWMD+bxldO+e0HuWyCjjDzMYSFoI+NTJybwHGmtmbgYeBs6L6xwKY2Z6Et+Ypkl7XmZtKdrns\nQPhWVFgM7Jeizg5mdr+kbxJ+470I3Gxmt7XZX8dpm3YFM+4u6fVLoTYvALcBjwCHAbuTvALaZgT7\n7Td0Jr1hzmjR5WJmS4Al0ecVkuYDI8zs1li1e4D3RZ+fAjaPvB2bAy8TjN2OkCToadeG2+DbIukN\nwKeBUcBzwM8lfcjMLq+uOzAwsO7zhAkTmDBhQsrLOk568ii+2YYzriUsZflrgoifCmzaRPvxwHnA\ns8CrM+xXOzwGLMr+tBkMikZu5XEEAY/zMUI8KGZ2s6QPE4R9M+DTZras/avXJum2ngBGxvZHEizw\nRnV2jMomAHeZ2bMAkq4l/K5rKOiO0wnyGQeeJU8RolcgxCFs38I5hhNEfQbrDcxes1O0Vbgjm9M2\n856rgaQtgKuB081sRaz8HOBlM7si2j+e8GC3J4QU/VbS7Wb2WHs9qE2SoN8HjI7eRE8CHwCOq6pz\nHTAJuErSeGCZmS2V9Efg/0oaDrxEcOD9odZF2p204Tj1KL+QvwRMBx4EDgH2or0J4PsRrPRngNe0\n3bvcUsflMmNx2BohaRhwDXCZmU2LlZ8AHE74h6jwVuAXZrYGeEbS74C3EH56ZE5DQTez1ZImATcT\nHsFFZjZf0inR8almdqOkwyUtBJ4njL5gZrMlXUp4KawF7gcu7MRNOE415RdyI0QC3wLsCnyS4KJt\nl02B/QlW+jEZnC+n1FG+CaPCVmFKlTNFkoCLgHlmdm6sfCLwOeAgM3sp1mQB8A7gMkmbE34Cfbvt\n/tdBZmnd5B3qgGQwAJQ88ZDTNcov5s8Q3CsvESYDjWxcvWlWEqz0jwDbZnzudhnAzJJGeBsiyezM\nlHW/xqDrSTqQMHI8h/VjjGcTHtjGrA8TutvMPilpE8IL4M2En04/MrNvttP/RuRqvpS7Xpx2KbeY\nv0zQkvuBgwi/3DsxQ2YTwiSj6ZR26kjrUS53UtunNbpO/ZWEQY2u4NkWnVJRToPACDHi3yUEjP0b\nwdfdyemObyHEPzzZwWv0kJJmW8ydoJfbwnKcZvkbcAVwOyE+/H3All247sbAgQRfeglxQXecYlAO\nK30VIUTvB8Drgf/D4PC9brAPYQ5NQthHEfFForuHW+lOf7OQkHvlKeAUgqXcC3NxGPA2Smmlu4Xu\nOMWhmFb6c8DPgOsJ+fCOpffZD8cRomr+0uN+ZIwLendxK93pH9YQFpv4PrANYcr+rj3t0XqGAm8n\nRLyUCF8k2nGKRTGs9D8DU4FHgZMIc1CG9bRHG7IXYbm6RT3uR4a4hd593Ep3yssK4BeEGeQHEUKV\n85IQq5ohhD5OJ32+vpzjgu44xSN/Vvpa4F7ge4Tke6cCY0lOb9tr9iC8hP7U645kQ0ldLrl/B/ns\nUac8PEGYsj8U+Cj5m1bfiLiVvjP5fwEl0Ga2xbySe0FvRNwl46Lv1MNsco/ddy8SJgYtICQdfTPF\nFMTdCYtJL6TOTPfiUGjlq08hXC61/jNWl7m/3ckfBswmLMsLwb2yF8UUcwhyMYFS+NLd5ZIPXLid\nVui+lf4IIev0xsAHCSs1loExhARhfyQsNl1QCqd86SiEhQ7pVkF3sXfyw0aEiUIfojxiDoOt9LW9\n7Uo7eJRLOXDR71+6Os6y5aWE6JWai3QVnDcRpGNBrzvSOi26XCSNlDRd0kOSHpR0WlT+DUnzJT0g\n6VpJr4zKR0l6UdKsaPteJ2+rdIJeT7DjFr6LutMxtpwcNoDNryEI+gu97FEHEHAwhbbSN025bcgq\n4AwzG0tYfehUSWMIS0eNNbM3Aw8DZ8XaLDSzcdH2yc7cUKB0gl5NGleN0z901ErfsurcG+1M8Dnf\n1blr9ozRhPGBh3rdkdZo0UI3syVmNjv6vIKQqH6Emd1qZpW32z3Ajh2+g5qUUtDTWOIu8v1JR/7d\n41Z5NZv/AphJWG63TFSs9BkU0krPwIcuaRQhe1nVyqN8DLgxtr9T5G6ZES1h1zEK6PZPhwu2E6cT\n34crbTbHveIXjStt9DpC/PbvgMMy70NveQNhYeq5hNj6AtGm8knaArgaOD2y1Cvl5wAvm9kVUdGT\nwEgz+7ukvYFpksaa2fL2elCb0gp6GpJmoVZEwCctFZdOvdivDL+6ufIfRyeL+ua/gOd3Afan46sN\njZsMs7plzFSs9OsIL60CBW7XUb4Z98KM+xo3lTSMkITnMjObFis/ATgcOKRSZmYvExaDxczul/Qo\nwV91fzvdr9s3s95OEJBkMNCz69cS61oi4KJenyxffFmmeuiUVV5NoqADLB8ffXhXth2qMK7qmXVN\n1AF+DOwJ7N2Faw1gZm3NzJJkNjdl3T0YdD1JAi4BnjWzM2LlE4FvAgeZ2V9j5dsAfzezNZJ2JgTx\n725my9q5h3r0tYUOgwXE3TTNkeXzyuu54tQSc0hrpU+D53cCDgBekV2nqoW8JxxMyBy5J4WRlNa7\neQAhNeYcSbOisrOB8wijxLcGzefuKKLlIGCKpFWEwYZTOiXm4BY6kH4WoVvpgUbPqtlnVO9crTzr\nbgt5nFRW+kufg1V3AO9uv1OQLOazBtiK21jGO+l8uoGfECJ63tLh62RkoadMGqmdaft63aQgr9PO\nkpUQ9ENmyLxa0p0Qc9tjIHzYE66ac1TDuolW+vIpBFF9kGDktbm0XArLfDPm8V7u5SpG8AJj27te\nIgcTls/bi0LISs66KOk7Kao9Z2ZfbHget9CbI43PPS+inrVvOy1J10t7rjT97pRVvk7MI5IEHRpY\n6cvjfbyNMNHoiFa7tp46oj78mals+8x57MMqfr7yEY5ha2YyhKXsx4sdtaAvB3YB9uvgNTKy0J9M\nWXdEdyx0SX8G/h/hrV9LlAV8wcwaJtDJ2XuqeOTV7x7vVzu/HPJqRXdLyCscu+e05q305bX6+Fbg\nO8CBwNYt9rIxL25zMk8P3ZpNnvgsAjZhNU9zCC+yW0eut56DgSsJg6N5W0avivwp37lmdkmjCpJe\nlXSSUk4s6iTxSUuN0gz0iixnxrZ6njTpjtu5fqfcK/XEvGmWT6kj5hBWKdqXEOzQJvUiWSRAbLpm\nGf+8yW5swkuVA+1fsyEjCInIEuL+coANSbd1rT9m386iTv7eUwUgj1Z5mkyUaa30PFnSeYpCSmOl\n1xfyOOMJVvqzdGod0U1XPsJVr7uYF7Z6L5stu5ZNF13QpYwyEwgDpPsQgj7yyZqcKl8U2vgpYBTr\n9dnMLJWPzn3oHaRbk5aK4JMuAmks9HqCfpz2avJqdwB/Bd7XZLsaNBognfVi9OHr7V8nNT8jWOud\nmOWejQ/9pZSZGDbdvLtRLpLmAD8kjJ5XciqYmd2Rpn1O31PlpZfWbyMrvZ+FvBmqrfTmhbzCeOB/\ngGeA12TQsyrWCXkvmECYbPRPwCY97Ed9Vm6S9tfDyx3tRw1eMrPzWm3sPvQOUj0wmaVvOUs/uYt5\nQHMHmqrfuphDELq3EpJbtUm1L72mmH++/euk5rWEhaSrc1blhzVDhqTaesB3JA1I2l/S3pUtbWO3\n0LuA+6TLRXtCHmdfgpW+FNi2vVPNmkJ3RTuJCcCPCPdYO7F4L1mT37wzY4EPE0KG4mksD07T2H3o\nOSGr2G2nfer50pu14NNxF/AX4NgMzpVG0JN86Z9PUSctvyBMoEqlRSnJxof+lL0yVd3t9Vy3feiP\nAmOipF5N4xZ6znEhzwedEXMIfua7CVlWR3ToGmnJ2sI/CPgBYaLRZhmfuz3W5Ff65gKvIvxsa5pE\nH7qkiZIWSHpE0pl16pwXHX9A0rhY+VaSro7W2psnaXyt9k52sdtO+8TFW3MHOijmECbgHEgmvvRU\nlnUt0f58VXlWwr41Yf3RuzM6X3asYUiqrQe8Clgg6RZJv4y269I2bviakjQEOB94J/AEcK+k68xs\nfqzO4cAuZjZa0n7ABYQhfAgOwhvN7P2ShhKy4Tt1qPi3Xch7T2dFvJq9CQtgLKa7K5d1w+f+duBC\ngiTk579/jn3otXyvqf3iDX3okvYHJpvZxGj/CwBm9tVYne8D083sp9H+AsJvrZeAWWa2c8MOuA/d\ncYB7gQWE8bB6vAf4ZYpzZSXUWfnSrydE9Ryawbmy8aHPt9enqjtGfy5UtsUkl8sOwOOx/cVRWVKd\nHYGdgGckXSzpfkk/kJQvR5rj5IZxhIlGf6lx7D3RVlTeRlhXdUVSxa6xhqGptm4h6fos6iQJelpT\nv/oNZgR3zt7A98xsb8IquV9IeT7H6TOGEn7YTq8qrxbyNMKelWWdlaX/SsLiF3dmdL72adWHLmmk\npOmSHpL0oKTTovJvRGOFD0i6VtIro/JDJd0naU70d72Qn7fFfOY1N0jOgZz0CnoCGBnbH0mwwBvV\n2TEqE7DYzO6Nyq+mrqDHv8SjCMa94/QbbwZ+CzwGnNbjvmTN24DvEnK8NDMz9jFgUea9acOHvgo4\nw8xmRwtFz5R0K3ALcKaZrZX0VeAsgt49A/yLmS2RNBa4mdoDJUemuPbKpApJgn4fMFrSKEJc1QeA\n46rqXAdMAq6KoliWmdlSAEmPS9rVzB4mDKw+VPsyWcapOk5RGQJ8GphG+JFbz3Wbxpf+dbKxsLOK\nS18S/T0DOKaJdjsx2MBLldIkkdUtCrqZLSG6GTNbIWk+MMLMbo1Vu4coSY/ZoOWu5gHDJQ0zs1VV\n553RUoeqaCjoZrZa0iTCW2UIcJGZzZd0SnR8qpndKOlwSQsJbpUTY6f4FHC5pI2BR6uOOY6zAauB\nvwF/oLMLRXSL54BfETTwaOCNve1ORBb+8cjQHceGOQ4+RkgMX837gJnVYp4luZopGp8t6aF7Tv8y\nDLgKuJjGOcw7HfHSjmW+Bvg9wW++H2HZvSwWvcgmyuVO26fmsftnLGfWjOXr9i+e8lTN60XulhnA\nV8xsWqz8HGBvM3tfVf2xwP8Ch5rZY+30vxG5mS6Vl2XbHKf3rCT82L2LIITdpl0XyyLgBsJg6El0\nKud7O7xcJ1f77hNeze4T1vf34ilPbVBH0jDgGuCyKjE/ATgcOKSq/o7AtcCHk8Rc0hHA9Wa2tlG9\neuQi22ItMXeBd/qXjYBTgO/TONAsTThjs+LcjpivIORvuZYwLvYh8ijmEHzoabZqJAm4CJhnZufG\nyicCnwOONLOXYuVbEd5uZ5pZmimzHwAWSvq6pIbrh9YiNxa64zhxnicEVPyGEM5YizQul7S0I+Rr\nCfETM4C9gFPJax70Cm340A8AjgfmSJoVlZ0NnEdYounWoPncbWafJASMvAGYLKlipR5qZn+tdXIz\n+1AU8ngc8OPgkuZi4EozW16rTZxc+NAbz1Z1X7rTr2wJTAUuZ/CP6VaEvJ4vvV33ymKCAbox8G5C\nLvROko0P/Xo7JLki8C+6vSczRSVtQ5g2/GlCdMxo4LykxS/cQnec3PIPQnDZdNa7ZfNilb8A/JqQ\nruBQwsShwsyQz20uF0lHAicQBPxS4J/M7Ololv08wi+BuuRe0D1ZldO/iBAVN5Xgn25nyCsel96u\ne+UB4DZgN4J7ZXgb5+sNrcahd4H3At82s9/EC83sBUknJTXOvaA7Tn8zmjCZ5iFgjzbP1a57ZSnB\nvbKaMODZ6/ztrfNyDn38UUba11eLeQUzuy3pHLmIcknCI16c/kXAOwgDjmt61IeVhLmFlxBcKydR\nZDGHfOZDN7PVwJooMqYl3EJ3nNyzMyGX+IOEfC/dwgi/DG6J+nAqecpp3g45drk8D8yVdAthoALA\nzCxVcp/CCLr70p3+RYS47uuA3aErYvRX4EaC7/79wOu6cM3ukeMl6K4hBPJXQv9EEwtc5PauHMeJ\nsxNh5uUDhKzUnWIVIePjfYQMifvSnRdId8lrlAvwqviEJQBJn07buBA+9AruS3f6m4MJE41Wd+j8\nfySkuH0W+D/A/pRRzCGfPvSIj9YoOyFt49JZ6O6accrL6wlT6WcDb8nwvMuAmwhulvcQJjaWm7xZ\n6JKOAz4I7BQtZlFhS8IbNhWFE/R6gh233l3UnfJyMPAzwhT7dv/7rgbuJiQB25+Qp7xwktASK/MX\ntngX8BRh9Y//Zv0srX8Ac9KepBT/eu6KcfqHHYFtCWt0tpMv/U+EQc+tgZOBV7XftQKRNwvdzP4M\n/BkYH+VZ38XMbotmiA4HEvO4QMF86BUqAm42ua6Yu8g75eVgQp7xVtZJWE5YDfJ/CYuIfZB+E3PI\nrw9d0snAzwnTgyG8wafVbzGYwlroLthO/zIC2IEQibJ/yjZrgHsJg6r7AEdAnZzg/UCO49BPJYQW\n/R7AzB6WlDrjWSEt9LS46DvlZQLwO+DlFHUfBy4kRLGcSEj01b9iDiEOPc3WA1aa2brFoKN0AKnj\n0Est6I5TXrYjTPb5Q4M6zxNcKz8DDgQ+Qhhzc1p1uUgaKWm6pIckPSjptKj8mKhsjaS9Y/U3lnSx\npDmSZkuql9y+wh3RMnabSTqU4H5JnWKz9ILuVrpTXiYQgiNWVpWvJQyafo9giZ9KSOxVnPS2naYN\nH/oq4AwzGwuMB06VNAaYS1gFuzqx1ieAtWa2JyHP8DejVY/q8QXgmeh8pxBGrr+Y9r4K60NPg4cu\nOuXmtYSY8XuAt0dlTxEyIkJYWGf7HvQr/6xs0eVkZkuAJdHnFZLmAyPM7HaAGlo9hpDQHjN7RtIy\nwiSCe+ucfw3BP3ZhK/0rpaC7kDv9w0HAjwhZEO8mJPA6hBCnXvof4C2ThX88Ci8cR3ij1uMB4AhJ\nVxJ8ZPsQIldqCrqk9wBfAkaxXp/NzF6Rpk+lE3QXc6e/2IZghX+PkLjrVGCznvaoCNQLSVw841Ge\nmPFoYntJWxDiP083sxUNqv6IYKXfR4gzv4vGeZDPJbhuHjSztYkdqaI0gu5C7vQvbydM29+n1x0p\nDPUEffsJu7L9hF3X7f9hyoZrSkgaRsiKeJmZNYwRj1won4m1/R3wcIMmi4GHWhFzKJGgO07/8vpo\nc9LSahx6NKB5ETCvOitivFqs/nBgIzN7PopaWWVmCxpc4kzgJknTWR+Tamb2rTT9K42ge/4Wx3HS\n0oYP/QDCaPMcSbOisrOBTYDvEHxgN0iaZWbvIuRp+JWktQTr+8MJ5/8yYTrvprQwWaA0gu44jpOW\nVqf1m9md1B9t3sD9YmaLgDc1cYntzezQFroGlGwY3GPOHcdJw8tsnGrrATdK+udWG7uF7jhO35Hj\nXC6fBP5d0susz76WOmyxVBY6uJXuOE4yec3lYmZbmNlGZrapmW0ZbanEHEoo6I7jOEnkOH3utZLe\nLaklbS5J1T4bAAAWK0lEQVSloLuV7jhOI/Iq6MAFwIeAhZK+KumNzTQupaA7juM0YjVDUm3dxsxu\nNbMPAnsDi4DbJd0l6cRoQlNDSivobqU7jlOPvPrQASS9GjgBOAm4HziPMA341qS2HuXiOE7f0aOQ\nxEQk/YIQt/4T4D1m9lR06CpJM5Pal1rQffao4zi1yHHY4g/M7MZ4gaRNzGylmSUm60l0uUiaKGmB\npEcknVmnznnR8Qckjas6NkTSLEmpV91wHMfpJDl2ufxHjbK70zZu2GNJQ4DzCcuDPwHcK+k6M5sf\nq3M4sIuZjZa0H2GUdnzsNKcD84At03YqS9xKdxynmh5FsNRF0vaE1b+HR0vYibCW6CtoIh9y0ito\nX2BhlI8ASVcBRwLzY3WOAC4BMLN7JG0laVszWyppR+BwwlvnMziO4+SAvAk6cBhhIHQH4Jux8uWE\n5F+pSBL0HQhLhldYDOyXos4OwFLg28DnCG+ZnuFWuuM4cfIm6GZ2CXCJpPeb2dWtnidJ0C3leaoX\n0pOkfwGeNrNZkiY03TPHcZwOsZJNet2Felwv6UOEJeiGELlezOxLaRonCfoTwMjY/kiCBd6ozo5R\n2fsIa+kdTsjt+wpJl5rZR6ovMjAwsO7zhAkTmDBhQpq+N4Vb6Y5TRB4jzK/JllYtdEkjgUsJK3Qb\ncKGZnSfpGGCAEHL4T2Z2f6zNnsBUwjji2uj4yjqX+F9gGTATeKnp/pnVN8IlDQX+SFh19kngD8Bx\nNQZFJ5nZ4ZLGA+ea2fiq8xwE/LuZvafGNaxRH7LGRd1xiswAZlbtEWgKSfYGezBV3Ue1+6DrSdoO\n2M7MZkfris4EjiKI+1qCcH+2IuiRhs4EjjezuZJeBTxXb4k5SQ+a2e6t3lvDsEUzWw1MAm4mRKr8\n1MzmSzpF0ilRnRuBP0laGN3MJ+udrtVO5o35dnFm57rTPJrTcbpNq1P/zWyJmc2OPq8gBIiMMLMF\nZlZrrdDDgDlmNjdq8/eE9ULviiz6lmhooXeDblvo4ZqtWenVQj5GJ7bch7iQH6gNfrg4jlOTbCz0\nHe2RVHUXa3Td60kaBdwBjI3EnWgt0LiFfjph2v5rgdcAV5nZN2qca270cQgwmuBrqrhlzMxSiXzp\nZopKUzqSxyUrq7yWRX6n/dJF3XG6SD0f+soZv+flGb9PbB+5W64GTq+IeR2GAQcCbwFeJCTbmmlm\nv66ql4kAlMZCr7a6k0Q9rZWeJOTNWOmN3Csu6I6Thmws9Nfan1PVfVqv3+B6UdbD64GbzOzcqmPV\nFvoHgHeZ2QnR/heBl8zsvxP6+FpCMAkAZvaXNP0tRbbFTg10ZmmVJ/nK3ZfuON1j5csbp9qqkSTg\nImBetZjHq8U+3wzsIWl4NEB6EPBQvX5JOkLSIwSXyx2EEJ+b0t5XoQVdmlJXzJNEPo1bJo31nST6\nAwxwW3KSNMdxusia1UNTbTU4ADgeODjKUTVL0rskHSXpcULakxsk3QRgZsuAbwH3ArOAmWbWSKC/\nAuwPPGxmOxEiDO9Je1+F9KEXIfRwgIGm27gv3XG6w5rVrcWhm9md1DeEp9VpczlwecpLrDKzv0ra\nSNIQM5su6X/S9q9wFnozYp5UN43oNmulD0R/qnEr3XHyw5rVQ1JtPeDvkrYEfgtcLuk8oNGg6yAK\nY6FnaZVPUVtjKjWZbxfzUzUeaLlNM3lnQkrjii/dLXXH6RyrV+Url0uMowjRMGcQ1hZ9BZBa/HIf\n5ZKFkFf85fWEPI2l3shXniTkcZIE3YXccRqRTZQLT6ScVb/Dpm1fr1UkvQZ4NmEi0iBybaHn3Spv\nRsgr1LPSXcgdp4v0xp1SF0n7A/8F/I0wMHopsA2wkaSPJgykriOXPvRG0Sutnq8RrfjSWxHzeriY\nO06XeWlouq17nA/8J3Al8GvgJDPbDng7QehTkTsLvRMRLK1EnNRijE7M5FwVK92F3HF6xOped2AD\nhpjZLQCSvmRmvwcwswWSUvvFcyPonRBymz+w7vOUMY3r1otOiR/PismTQC7mjtM78ifocdFuOm1u\nhVy4XDot5u2StZgD2KTszuk4TpOsTrl1jz0lLZe0nDCzdHl8P+1JcmOhZ0U9IZ88v3krvRNC7jhO\nDljV6w4MxswyGaUtlaBnZZV3S8ht0gA6P7trOY6TkjW97kBnyIXLJQvSiPnk+YlV3Cp3nH4gfy6X\nTCiNoGvMQK+70BLuS3ecHvBSyq1glEbQ05LGSs+KKed371qO4zSBW+j5x610x3FS4YLeW6Zn6Nvu\nJm6lO04OaVHQJY2UNF3SQ5IelHRaVH5MVLZG0t6x+vvG8qbPiVYw6hiFEPSKmKcR9SQrPSl0sROk\nEXWbPdDxfjiOE7Eq5Va75RlmNpawmMWpksYAc4Gjgd9U1Z8L7GNm44DDgO9K6lgimVyHLWZplfdC\nyFNxUq874Dh9SIthi2a2BFgSfV4haT4wwsxuB1BVEkAzezG2Oxx4zsw6FjSZSwt9OgN1xbxZK33K\nmHyI+QZW+klsIOZupTtOl8jAhy5pFDCOhCXiIrfLQ4S1RD/TVr8TyJ2FnpVVrjGNc7P0FLfKHae3\n1AtJ/OMMeHhGYnNJWwBXA6ebWcMVhczsD8BYSW8CfiVphpk911yH05GLBS6omnJ/UIp2B3cpkVbW\nTJ6dXEd7DXS8H45TTDJa4OK7KXXvVG1wPUnDgOuBm8zs3Kpj04HPmtn9da59O/B5M+vImpS5sNCz\nEuA8C7njODmixZBEBSf5RcC8ajGPV4vVHwUsNrPVkl4PjAYeae3qyeTSh35HijrVrpmiiPmUvZLr\nuC/dcTpM6z70A4DjgYNj4YjvknSUpMcJkS83SKqsMPQ2YLakWcDPgZPN7B+duq1cWOitMp2BVOLv\nOI4ziBazLZrZndQ3hKfVqP8T4CetXa15cmmhQzorvahi7la64/SYNSm3glFIC72oQt4MPjDqOB2k\ngIm30pBrQb+DwREvZRLyKXvVjnhxIXecLlDAPC1pyLWgxymTmNfChdxxukjOVizKitz60CvcQXnF\nvOJLdzF3+g2zyb3tQEl96LkX9LwxeY+wZcEAAy7mTt9REfOeinpJ0+cWxuWSB7IScihO3LzjZEXP\nrfI4BRTrNLigp6CWkE/eA6bMbf5cLuROv9FIyM0mI03pYm8i+tmHLmmipAWSHpF0Zp0650XHH5A0\nLiqrmQy+miwt3yzJ2r3iYu70G2ms8p5Y7itTbgUjUdCjZOznAxOB3YDjooTu8TqHA7uY2WjgZOCC\n6FC9ZPC5Z7ZdmVgnrdi7kDv9htnkfLlYquljH/q+wEIzWwQg6SrgSCC+3PIRwCUAZnaPpK0kbVsv\nGXxVW6B1F0bWpBHytLiQO/1IK0LedddLSV0uaQR9B+Dx2P5iYL8UdXYEllYK0iaD7xW1hPzoOVfy\niz2Pa9iu1ovIhdzpR3JtkVdTwJDENKQR9LQJ06tzFK9rl5QMfqAi+6+FRU/DTikvmBVulTv9RpYW\ncVZCXrtPjwGLMjn/IAroTklDmkHRJ4CRsf2RBAu8UZ0do7JKMvhrgMvMbINsZAAD267fksR88qQU\nPU7JwFRLFPOj56Tzpfugp1MUsowD77xVvhNwcGzLiJL60NMI+n3AaEmjJG0MfAC4rqrOdcBHACSN\nB5aZ2dKUyeAH0WigsSLm7Yr6wFRjYGo2KzVp7gCaO5DJuRynk2Q5UNmpQc+uuW1WpdyqqBe5J+mY\nqGyNpL2r2pwVRQAukHRYJ28r0eUSrbQxCbgZGAJcZGbzJZ0SHZ9qZjdKOlzSQuB54MSoeSUZ/Jwo\nwTvAWWb2q2Y6mbVVHmfahcdy1MlXNWxTz5fuQu4UgSzjwAvlJ29E6yGJlci92ZEreaakW4G5wNHA\n1HhlSbsRjODdCGONt0na1czWttyDBqSaWGRmNwE3VZVNrdrfQHYTksHXpTLQ2EjIJ0+CKeenP2dW\nFjm4kDvFIUsB7paYdyXipUV3Sr3IPTO7HSA4JQZxJHClma0CFkVG777A71vrQWNyO1O0k1Z5NWmt\n9Gn6Y3adcpwO0oz4JgloaazyOBmELaaM3BvBYPFeTLDUO0LuBH2d9TsXbNJAw7pJVnqWg5Qu5k4R\nKKJFXu/aHbXS64UtLp8BK2YkNk+K3EsgO3dBFbkR9CzdGBsI+SlTYGrjL2c9K92F3CkK7QhwtYCW\n0iqPU8/lMnxC2Cos2fClkiZyL0bdCMBOkIv0ufXEXOfXLo9T7Zpxq9zpR/IevZI7WgxbTBm5F3ek\nXwccK2ljSTsBo4E/ZHMTNS5s1jHrP10HJKOBCCe5XdadJ4X4J1npAEedfJULuVNIyiTE9d0tA5jZ\nBiOPzZ1bxi4pdW+hBl1P0oHAb4A5rHednA1sAnwH2AZ4DphlZu+K2pwNfIzwijjdzG5up/+NyL2g\nQ2NRTyXkFZIE/ZQepPF0nIwog6An+80zEvSRKXXvcbV9vW6SC5dLqzQl5tBYsF3MnYLTk7ziGSFN\n6W7/SzpTNDeDoo3Q+QODrPSmhbwRLuSO01N8gYvsKISgV8hEyCsRLy7kTgmRppTC9dJxSpptMSc+\n9MuAhT3th+OUhSIKenorPSMf+pYpdW+5+9Adx+khRfald42S+tBzJOi79LoDjuP0iK7/qmgx22Le\nKZQP3XGcZIrocuk6BbS+05AjC91xnHYpspgXue95IUcWug+KOk6ruBg6kBsL3cXccVqhbLlXynQv\nvSBHFrrjOM3g4tcOBRzxTIELuuMUjLILeVdWLCrpqKgLuuMUhLILeXcpp4WeEx+64ziN6Dcx7/z9\nvphyKxZuoTtOjuk3Ie8erVnokkYClwKvJeRDv9DMzpO0NfBT4PXAIuBfzWyZpA8B/x47xZ7AODOb\n00bn6+IWuuPklH4X887ef8tz/1cBZ5jZWGA8cKqkMcAXgFvNbFfg9mgfM7vczMaZ2Tjgw8CfOiXm\n4Ba64+SOfhfy7tCahW5mS4Al0ecVkuYDOwBHAAdF1S4BZhCJeowPAhsuXJwhLuiOkxNcyDekcxEv\n7Ue5SBoFjAPuAbY1s6XRoaXAtjWa/CtB+DuGC7rj5AAX825Tz0K/F7gvsbWkLYBrCGuELg9rRwfM\nzEJa8EH19wNeMLN5LXc5BS7ojtNDXMh7Rb0Ilt2jrcL3N6ghaRhBzH9iZtOi4qWStjOzJZK2B56u\nanYscEV7fU7GBd1xeoALeTo6N8GoNZeLgil+ETDPzM6NHboO+CjwtejvabE2GwHHAAe22tu0uKA7\nTpdxMU+m8zNFW55YdABwPDBH0qyo7Czgq8DPJH2cKGwx1ubtwF/MbFGrF02LC7rjdAkX8mS6t9pS\naxa6md1J/XDvd9ZpMwN4a0sXbBIXdMfpMC7k6eju0nnlnPrvgu44HcTFPJnerIHqybkcx0mJC3ky\nvV3M2i10x3EScCFPR2/FHIqYeCsNLuiOkxPiIlfWF0PvhbyCW+iO43SAWiInTSmtqOeDcvrQE7Mt\nSpooaYGkRySdWafOedHxBySNa6at45SJZi3Q/Fis/caqlFuxaCjokoYA5wMTgd2A46JUkfE6hwO7\nmNlo4GTggrRtnWoe63UHckS5n4U0JVHMK8dnzJjRhR51h/z86mg5fW6uSbLQ9wUWmtkiM1tFSP14\nZFWdIwjpIjGze4CtJG2Xsq0ziEW97kCOWNTrDrRMWqFOS5kEHfIi6n1ooRPy/D4e218claWpMyJF\nW8fpG9JY5bXaOJ2gPy10SzheQclVHKd/iAtxK0JednpvpZdzTVGZ1ddsSeOBATObGO2fBaw1s6/F\n6nwfmGFmV0X7Cwgrd+yU1DYqT/vScBzHwczaMiCb1Zx2r9dNksIW7wNGRytzPAl8ADiuqs51wCTg\nqugFsMzMlkp6NkXbQj0sx3GKT5k1p6Ggm9lqSZOAm4EhwEVmNl/SKdHxqWZ2o6TDJS0EngdObNS2\nkzfjOI7TzzR0uTiO4zjFIXFiURZI2lrSrZIelnSLpK1q1NlU0j2SZkuaJ+m/mmlfFFI+i5GSpkt6\nSNKDkk6LHRuQtFjSrGib2N07yI4MnkVffS+iej+StFTS3KryvvpeRPXqPYvSfC+apSuCDnwBuNXM\ndgVuj/YHYWYvAQeb2V7AnsDBkg5I275ApLmXVcAZZjYWGA+cKulN0TEDvmVm46LtV13pdWdo91n0\n2/cC4GLCZL1q+u17AfWfRZm+F81hZh3fgAXAttHn7YAFCfU3Iyy/vVsr7fO8tXIvhPUJD4k+TwY+\n2+v7yMmz6MvvBTAKmFtV1pffizrPojTfi2a3blno25rZ0ujzUmDbWpUkbSRpdlRnupnNa6Z9QWjq\nXqIooXHAPbHiT0V5cy4q+M/Jdp9F334v6tCX34sOtC8smWVblHQr4W1YzTnxHTOzenGgZrYW2EvS\nK4GbJU2wsB5fqvZ5IYtnEZ1nC+Bq4HQzWxEVXwB8Kfr8ZeCbwMfb7nSH6PCzSN0+D2T1LOrQl9+L\nJIrwvciSzATdzA6tdywauNjOzJZI2h54OuFcz0m6AdgHmAE01b7XZPEsJA0DrgEuM7NpsXM/Havz\nQ+CX2fU8ezr5LOjD70WDc/fd96IBhfpeZEm3XC7XAR+NPn+U4AcdhKRtKj8TJQ0HDgVmp21fINI8\nCwEXAfPM7NyqY9vHdo8GBo3wF4y2nkWa9gWirXvpt+9Fh9sXl2446oGtgduAh4FbgK2i8hHADdHn\nPYH7CSI+B/hcUvsibimfxYHA2uhZzIq2idGxS6Pn8wDhi7ptr++ph8+ir74X0f6VhJnXKwnJ707s\nx+9FwrMozfei2c0nFjmO45SEbrlcHMdxnA7jgu44jlMSXNAdx3FKggu64zhOSXBBdxxnHZKOiRKh\nrZG0d506jRKmfUPS/GjG6rXRJEEkbSzpYklzFBLwHZSiL5dLWiBpbjT7NbN5M2XFBd1xnDhzCXHs\nv2lQp1bCtDHRsVuAsWb2ZkLY4FlR+ScIK5btSZhj8s1ojkEjLjOzN5nZHsBw4KSW7qiPcEF3HGcd\nZrbAzB5OqLPEzGZHn1cA8wkx4pjZrRZSeEDIubNj9HkMMD2q8wywDHgLgKTDJN0laaakn0naPKp3\nU+yy98bO5dTBBd1xnJapkzyuwseAG6PPDwBHSBoiaSdCWo8dJW1DyN9yiJntA8wEPlN1jWHA8UBc\n4J0auE/KcfqMBomxzjaz1DlgGiVMk3QO8LKZXREV/Yhgpd8H/Bm4izADeDywG3BX5IHZODoW53vA\nHWb2u7R961dc0B2nz7AGibHS0iBhGpJOAA4HDoldcw0xy1vS74A/ArsQFqP4YJ3rTAZebWafaLfP\n/YC7XBzHqUfNQcuE5HETgc8BR1pYhaxSPrziG5d0KLDKzBYQXDUHSHpDdGxzSaOjzycBhwE1xd7Z\nEM/l4jjOOiQdDZwHbAM8B8wys3dJGgH8wMzeLelAQhTMHMLSdwBnmdmvJD1CcJv8LSq/28w+Gfna\nf0VwsywGPm5mj0fXPBj4GrBJ1OYcM7te0ipgEVBx51xjZl/p0K2XAhd0x3GckuAuF8dxnJLggu44\njlMSXNAdx3FKggu64zhOSXBBdxzHKQku6I7jOCXBBd1xHKckuKA7juOUhP8PF7E+7mg0AmEAAAAA\nSUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f2bf9fa0410>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.pcolormesh(lons, lats,bathy,vmin=205,vmax=234)\n",
    "cbar = plt.colorbar()\n",
    "cbar.set_label('Bahtymetry [m]')\n",
    "plt.plot([lons[j-stp,i-stp], lons[j-stp,i+stp]], [lats[j-stp,i-stp], lats[j-stp,i+stp] ] ,'k')\n",
    "plt.plot([lons[j-stp,i-stp], lons[j+stp,i-stp]], [lats[j-stp,i-stp], lats[j+stp,i-stp] ] ,'k')\n",
    "plt.plot([lons[j+stp,i+stp], lons[j-stp,i+stp]], [lats[j+stp,i+stp], lats[j-stp,i+stp] ],'k')\n",
    "plt.plot([lons[j+stp,i+stp], lons[j+stp,i-stp]], [lats[j+stp,i+stp], lats[j+stp,i-stp] ],'k')\n",
    "plt.axis([-123.2,-123, 48.7,48.8])\n",
    "plt.plot(lon, lat,'r*')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No where in this neighourhood do we get close to 234m depth. Is the depth record in the ONC file accurate?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# IOS Cast"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('> open ', 'strat/comparisons/SOG_2000.dat')\n",
      "> load CTD_DAT\n",
      "> load STANDARD_KEYS\n",
      "('> close ', 'strat/comparisons/SOG_2000.dat')\n",
      "> complete\n"
     ]
    }
   ],
   "source": [
    "ACTDR.load_dat('strat/comparisons/SOG_2000.dat')\n",
    "data = pd.DataFrame(ACTDR.CTD_DAT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "lon_min=-123.49\n",
    "lon_max=-122\n",
    "lat_min=48.75\n",
    "lat_max=48.8\n",
    "data_BP = comparisons.isolate_region(data,lon_min, lon_max, lat_min, lat_max)\n",
    "data_fall_BP  = data_BP[(data_BP.Year>=2014) & (data_BP.Month >=10)  ]\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(-123.027, 48.771999999999998, 211.08877605139131)\n"
     ]
    }
   ],
   "source": [
    "lon = data_fall_BP.Longitude\n",
    "lat = data_fall_BP.Latitude\n",
    "depths = data_fall_BP.Depth\n",
    "print(lon.values[0], lat.values[0], max(depths.values[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(343, 292, 131.0)\n"
     ]
    }
   ],
   "source": [
    "j,i = tidetools.find_closest_model_point(lon.values[0],lat.values[0],lons, lats, bathy)\n",
    "print(j,i, bathy[j,i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model is 80m shallower! Neighbouring points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[186.0 101.625 62.375 63.3125 59.0 128.0]\n",
      " [208.0 178.0 116.625 55.375 45.4375 94.9375]\n",
      " [204.0 189.0 159.0 120.0 102.625 114.0]\n",
      " [205.0 180.0 201.0 131.0 129.0 149.0]\n",
      " [192.0 200.0 197.0 140.0 142.0 172.0]\n",
      " [139.25 164.625 192.0 188.0 174.0 184.0]]\n"
     ]
    }
   ],
   "source": [
    "stp = 3\n",
    "print(bathy[j-stp:j+stp,i-stp:i+stp])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f2bfc152610>]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAETCAYAAADETubIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXe8FNXZgJ9XbKixx4ZEbKBGjRWxfWAHC5YYI0aN3cTP\nhhoVk8jFfPlsUbF9sWAjCoqxoQFEVOwNRMUIKooKGkFREUWQ8n5/nFnu3Lmzu2d3zszO7j3P7zc/\ndmbPnHmHvffZc985RVQVj8fj8dQ3S9U6AI/H4/Ekx8vc4/F4GgAvc4/H42kAvMw9Ho+nAfAy93g8\nngbAy9zj8XgagLqQuYicKyKLRWT1Iu/3E5F/i8hEERkiIssFx+8VkQnBNlVEJlhca5SIfC0ij7q+\nD4/H40mL3MhcRHqIyB0xxzsC+wAfFzmvE3AysJ2qbgW0A44EUNUjVXVbVd0WeCDYynEFcEw19+Dx\neDy1IjcyB4qNXroaOL/Eed8CC4AVRGRpYAXg03ABERHgCGBosN9ORK4UkVdF5E0ROWVJEKpPAd9V\nfxsej8eTPXmSubQ6IHIwMF1V3yp2kqp+BVwFfAJ8BnyjqmMixXYHZqjqB8H+iUG5rkBX4OSghe/x\neDx1ydK1DkBEXgaWA1YCVg/ltfsD/YB9w8Vjzt8YOBvoBMwG7heR36jqPaFifYAhof19ga1E5PBg\nf2VgE+CjpPfj8Xg8taDmMlfVbgAi0h04TlWPD/a3BDYE3jRZEtYHxotIV1WdGapiB+BFVZ0VnPcg\nsAtwT7C/NHAosF3k0qer6hPFwnJxbx6Px5MVZdMsItJTRCaLyPsickHM+5uJyEsiMk9Ezq3k3GhV\n4R1VfVtV11bVDVV1Q2A65iHnzMh5k4FuItI+yI3vDbwTen9vYJKqfhY69jhwWiB6RKSziKxQLBaP\nx+PJOyVlLiLtgBuAnsAWQB8R2TxSbBZwBvC3Ks4No5RuES95T0TWE5F/Aajqm8BgYBxQyK3fEjrv\n1wQPPkMMwgj/dRGZCPwd0wsGEXkOGAbsJSLTRGSfEjF5PB5PLpBSU+CKyM5Af1XtGexfCKCql8WU\n7Q98p6pXVXqux+PxeJJRLs3SAZgW2p8eHLMhybkej8eTO0Sko4g8HQxSfFtEzgyO/yo4tkhEtouc\n0y9INU8WkX3ja05OuQegSR4E+oeIHo+n0VgA9FXVN0RkJUynjCeAiZiOFjeHC4vIFphU7xaYxuwY\nEemsqotdB1ZO5p8CHUP7HTEtbBuszhURL32Px2ONqibqoFCpc8LXU9XPgc+D19+JyCRgPVV9Mqg7\nevrBwFBVXQB8JCJTMGNbXq7+DuIpJ/NxwKbBgJrPMN8wfYqUjd6F9bm1XrpOZEAKtT6N6thENUTj\nUu3vpJ7qeRrYw1FdecLfF0D72WeWfP+HVa5LGE916INNS17LYW7q/B/Lcn8q8V7gtm2BV0oUW4+W\n4k4t3VxS5qq6UEROx3TlawfcpqqTROTU4P2bRWQd4DXMwJvFInIWsEXwrdXq3DRuohrSEXjra1Qq\n4GJxVSvywrlZ3K+n/ign8GjZrIQeFngaLFPk+AfAhxbnBymWfwJnqWql03+k0notO2hIVUcCIyPH\nbg69/pyW6ZSS59aKWsmsnNBt46rmiyFMUqGr9qepSWlq6r8kHk/9UYm8i52fhtDTlneUYuLrEmwF\novOCAIjIMphJ++5W1YfLXCqabl6fyNxRrqj5CNA0qa1wOi15FRVxtXFlLfTotXr06FF1XfmmU60D\nSIlO5p8bw5/j14lrLXwhJJV6E00A9H8waUSV077K84KBibcB76jqwGLFQq+HA0NE5GpMemVT4NUq\nL186ttrnq0VdxtA4gilOEqFD6f8jV2khT425sfjn2P7o5EIvUKnQCwKPYit0OczNA9BbLcueTMvr\nichuwLOYAYoFcV2EmV/qemBNzBxRE1S1V3DORcAJwEJMWubxJPEXo+5l3lZl4lLoaX45eDKihLzj\nyEroxeQdh43QXcm81cIJRTie5NfLirqUuZeHIamEXeI/kxpQocCjpCX0SgQepZzQXcn8bsuyR+Nl\nbh+Apcy9LOLxQm9jHBR83j3dVOdK6J+vvTrXzHNSVUmhu5L5vZZlj8TL3D6AEjL3crDDC73BOajI\n51tjoX++dusledMWuiuZ26wfCfBLvMztA/Ayd4IXegNRTN5xZCj0OHnHkabQXcl8uGXZ3niZ2wdQ\nJs3ixWCPF3odU4nAo6QodFuBR0lL6K5kXmxVmij74GVuH4BFztyLwR4v9HqiJxy0k7OqXDH7XDfu\nci30AYdBE25k/rRl2T0cXC8r6kLmppwXgy1e6HmliHFrLHT9sKWrvr3RQSwBroReoAk3Mn/esuxu\nDq6XFWWXjcsLeRJU3smTQP3n1jO0FeHRUvM0VcAo+6L6oSzZoqz8327CAei7vLu6XLK05VZP1E3L\nvLl8fkSVd/Ik0rb1uVXZRE65hR4n7lLksYXehJuW+XjLsts7uF5W1J3MzTltSQzJ8ELPCkdJa4dC\nr1TeceRN6E24kfmblmV/4eB6WVGXMjfnNbIY3OKFnhW1F/q8e5u9s9zFLoIxJJX6ysFiaQNs+wQW\noQk3Mn/HsuwWDq6XFXUrc3NuI4vBLV7oWZG90MMCj1JLoa9cZLXLJEJvwo3Mp1qW3dDB9bKirmVu\nzm9kMbjFCz0r0hd6KYFHyVLoxQQepVqhN+FG5rZrX67v4HpZUfcyN3U0shiSkSeBR2nsz82t0CuR\ndxxpCd1W3nFUI/Qm3Mj8c8uy6zi4XlY0hMxNPY0shsrIs8CjNPbnllDoAwKRn+fGJc6Evg1wv5uq\nKhV6E25kPsuy3+EaC73M7QNwuDhFY4uhNPUk8CiN/blVKPQB8amVmgt9m5hjNRB6E25kPtuy//sq\n87zM7QPwKw1VTT0LPEpjf24lhF5E3nFkKvQ4eceRsdCbcCPzuSvalV3hey9z+wAcy9zU2ZhiaCR5\nx9Gon5shJPQKBB7FhdDnzYMXT4Y/Pg+X7AIHbhy8YSvwKBkKvQk3Mv9xFbuyy872MrcPIAWZm3ob\nQwyNLvAojfK5RWk/+0x+GLiak7qqEfoXX8DIUfCvEfD0WNhyS+i4Pjw1EkZfANtskDCojITehBuZ\n61qWZWd6mdsHkJLMTd31KYa2JvAo9fq5RSmsYh8mK6GrwnvvwWP/gsdGwL//DXvtCYfsBfvvDWuu\nYcrdfwmcPhhGnAfbb5gwqARCHxSS+KclyjXhSObrWZb9rNWCzh2BwcBamAWdb1HV60RkdeA+YAPg\nI+AIVf0mOKcfZkHnRcCZqjo6SfxFY21kmZv68y+Gti7vOOrhc4sSJ+840hL6woXw4kvwr0Dg8+bB\ngQfAIXtCj11g+SIP/R76K/zuDnj0HOi6cXwZayoQ+qASrfBiQm/Ckcx/Zln2k1YyXwdYR1XfEJGV\ngPHAIZi1n79U1StE5AJgNVW9UES2AIYAOwIdgDFAZ1VdnOQeYmNtdJmba+RPDF7g5cnj5xbFVuBR\nXAl95inC6CdMC3z0E7DBz+CAA+CwPWGbLUEstffYZXDCrfBIX9h504RBlRB6KYFHiRN6E45kvpFl\n2Q9LX09EHgZuCLbuqjojEP5YVd0saJUvVtXLg/KjgCZVfTnJPcTG0hZk3nyt2srBC7xyav2ZxVGt\nwKNULfRvPoZ3H4V3h8P0l9mv+xzTAu8O61umD+IYdSUcexM8cBbs3qX6eoAWQq9E4FGiQm/Ckcwt\nv7Dk/eLXE5FOwDPAlsAnqrpacFyAr1R1NRG5HnhZVe8J3hsEjFS1XobUmnqbsjcRqv0zl4MXeDJq\n8ZmV44dVrnMi9PZnf20n9MWL4T/jYfJwI/A5n0HnA/jdFT+yxb4HsPxKy3DyF/9IHE/PP8A9Aodd\nC/efAT02r66e+YcBh8E/LLv/laIDpXPoVVPEfGO/h7Fzy58epFgeAM5S1TkS+hNIVVVESrVQU2m9\ntqmWefM105ODl3c65E3okHILfcEPMPUpI/D3HoXlVoEuvTn/8jfZqNuaLNWu9boyLoQO8PTzcMRv\nYehpsPeWdufMPyz+uAuhQ7PQm3DUMrfshilvtL6eiCwDPIZpYQ8Mjk0Geqjq5yKyLvB0kGa5EBPz\nZUG5UUB/VXW0IkkorrYoc3Ndd3LwAs+Ghhf6dzPgvX+Z1vfUp2HdbaFLb/5y/bOsvenKVvW4Evqz\nL8HhR8Pg30HPrePLFBN4FJdCb8KRzHewLDuu1QNQAe4CZqlq39DxK4JjlwcCXzXyALQrzQ9AN0lD\nem1W5uba1cvBC7w2NJLQVRWdPJlFI0axaORIdOL7sPG+nHDhZ2zZqwMrrbFcVfW6EvqLr8IhfeD2\nk+HAbe3lHYcroZ+MI5l3syz7ciuZ7wY8C7xFc7qkH/AqMAz4Ga27Jl6E6Zq4EJOWeTxJ/EVjbcsy\nN9e3k4OXd36oZ6HrggUsfullFo0cyeIRo2DRQpbq1YuRnW6l+9Zw1wHHOInHldBP+/DX/PPA+7jx\nBuh9ULK6XAjdmcx3tSz7gh80ZB9AjWVuYoiXgxd4fqknoes3s1k0ZgyLR45i0ZgnkY02ol2v/RjX\n4VK23qh198Fb962t0C/76dkt9qePm8H9+w9l4DVw2KHJYkoi9Kv1dSbJdm5k3t2y7DNe5vYB5EDm\nJg4jhzwIPCyqPMRTIE//R5BvoS/+6GMWjxrFohGjWPz6eJbaZRfa9erJR2udw3prlK8na6FHBR7l\nszdmMmy/IfztSvjV4clishX61fp6q2POZL6XZdknvcztA8iJzGtNKTHVUp55/aslfyJfDHwGvIts\n8SX6xZe0229fHtj4HvbeDla0nHI1TJpCLyfvOD6f+CX37XM3//tXOKpPspiKCT1O4GGcyXw/y7KP\ne5nbB9CGZV6JkLKUZ16fI+RP4D8CHwLvAu8BKwBdePG05+l6ILRrl/wKroQOMAuLPwnKMOOdWQzb\n+x809YdjE4ZWEHo5gYdxJvPelmWHe5nbB9DGZJ7XHjTVxpXHmNJlDkbc72I6LawHdGHK+Y+zcdSV\n1U4pGyGJ0E9ZefCS15d+27dESXtmTv6KYXsN5o9/hBOOq66OXVd4HoC5rFDRec5kbtkzRx5sIJmL\nSE9gINAOGFSYYyBS5jqgFzAXOE5VJwTH+wFHY/4GnQgcr6rzI+c2tMxdC8mVPPPazz5/AldgJkbe\n7wKzgE2458i36dUFVivnooyFHpZ3HK6E/uX7X3Pfnndx/nlwysnlyxfkHUclQncm819Zlr2/QWQu\nIu0wP8F7Y/rsvwb0UdVJoTL7A6er6v4ishNwrap2C+YteArYXFXni8h9wAhVvStyjYaTedpCqlae\neR35mj+BLwQ+oVngAF0Yc/Ir7N4Jlq10EoyUhV5O4FFcCf2rD2dz3x53cPaZcNrvW79fSuBRbIXu\nTOaWOX8Z2jgy3xkz9LRnsN9iaGpw7CbM0NX7gv3JQHdgAfAS0A3zt+lDGNGPiVyjIWSe1zlfsoyr\nEqHnT+A/AO9j5P0BsAbQhbfOfoot17GffbAojoVeqcCjuBL61x/NZtged/D7U+GsMysTeBQboTuT\nuWXmSv5RPzIv18boAEwL7U8HomtexZXpoKqvi8hVmCbOD8DjUZHXM7WWkciAWHnWMq5iMRXeyx9f\n0dz6/gzYEOjM9v8ZwLLrrA7AqfyJF6/YM/ml3iCR0GW/UIPnJ4mjod/K1zgR+hV/vR26DeDi23bh\nOjmQdS6ovq4VmFtxDr1qHDyczhvlZG7bZG71zSUiGwNnA52A2cD9IvKbwlSQYZqampa87tGjBz16\n9LC8bLbkTUgFeeYprrDQ8xSXYTEmW1gQ+FygC12Gn8Qqe21HuxXi+w/ucv5TNRF6C4GHmYMzoUPl\nrfR+J1/T8sBKHZm/48t8eccu6I8LWffPx1UdU1To348dx9yx46qurygNOF9suTRLN8xE6oU0S4uJ\n1oNjN2EmYr832C+kWXoA+6jqScHxY4BuqvrfkWvkOs2SPyF5KuNHTNqk0H1wJaALW758JCvt2AVZ\nqvXsg8VwInQoKfSiAo/DgdALlBN6K4HHMfdzlp+wM6v9ag/WbToBSZCbKtZCd5ZmOd2y7A31k2Yp\nJ/OlMb8Fe2H+Dn2V0g9AuwEDgweg2wB3Y5ZLmgfcCbyqqjdGrpE7mXuB1zvf0tx98GNMJrAL2354\nDMtvuG6imtMQekUCj5Ki0K0EHuWHmSz/5s6setCurPfXU5wL3ZnMz7Ise239yLzkHxuqulBETgce\nx2SZblPVSSJyavD+zao6QkT2F5EpwPeYtfBQs0beYGAc5u/b14FbUryXqvHyrncUmEFz+uQrYFPg\nz5iOWKuys6NVupykXF4wmwx30IhxlHIB6CerwkkJu5m2X4t527zC7BE7owsW0uGK06oWeiHlMqn3\ntsliiqOtpVkyCaAB5jP31IKFmEE7BYG3A7oA/YFdgGVaneFK6FBFC/2F+MNOhA7VC31OzO9BUqED\nzPuKFd7ZmZV2/wXrX3NmxUIfL6G+EgeF4nlU3LTMLR/UyuX10zJvkzL3Iq9XvsKkTz7B5MF/ihH4\nFcBmxDyHj8WV1MsKvYjAo2Qu9DiBR3Eh9Plfs8LkXVmx6+Z0vO7sss8nWgg8SkHormR+kWXZ//Uy\ntw/At8w91jyKWROgByZjt3bVNaUidEt5x5Gq0G3kHYcLof84mxXf3432W23Ez/5+Xguhl5R3HAf1\ndyfziy3LXuJlbh9Ana40lJhDY35RHqr9F8xR2rHF/hCZVqRktugbTQDINlcCNwAjMb1eq8eF0F/c\nek/TPd0BzoQOgKOfJSdCn8NKH+7Gcp07Muv27QH7HkStaXIj80ssy17sZW4fQJ2sNOSEOIFHqYHQ\nowKPUiuhFwQeRba5BrgaGAFsnOga1Qj9xa1j0iu5EHpKPztJhT5oAKaL6BBgZeAQqhe6I5m3mmGq\nSNkLvMztA8hB18RUhW4j8CgZCL2cwKNkJfRiAo8i21wPXIYR+qaJrmkj9FiBR6mJ0DP68q9U6IPi\n4voRuBczVfChVDcM05HMr7Ise66XuX0AOZC5icPRL0U18o7DsdArlXccaQjdVt5xyDY3AZcAj2Ee\ngFZPVOhW8o4jdaHXMBVXSuix8o5jAXAfsCzwSyoXuiOZX2tZ9qxWCzrfDhwAzFTVrYJjvwBuAlbE\ndLH6jarOCd7rh1nMeRFwpqqOThJ7yVhrLdK8yBwaT+j6fNOS17/Z9bZksQS4ELre0AS7JY8FQLYZ\nBPwJ83D054nq0q0czQviXOi1f5ayhLDQrQUeZQFmIft2wOFU1unbkcxvsCx7eiuZ7w58BwwOyfw1\n4BxVfU5Ejgc2VNWLRWQLTG5pR8zItTFAZ1VdnCT+orHWWqR5kjnUv9DDAo9SS6HrDU2tDzoT+p3A\nBcAjwNYVnauXRwR+t5uYEgv9DfOPfNKUNJKcshC4HzPg6wjshe5I5jdZlv1d6zRLML33oyGZf6Oq\nqwavOwKjVPXn0elPRGQUZnoUdwMewnHVWqR5kznUn9BLCTxKVkKPlXcczoR+N3AOZqbl4iMGW8k7\njloI/Y3ibzWu0BcB/8S01H9N3ECv1jiS+SDLsidZyfwF4ApVfUREzsEIe2URuR54uTC5oIgMAkaq\n6gNJ4i8aa61FmkeZQ76FXom840hL6NYCj+JM6EOBszCC2GHJcSuBR8lC6CUEHqWxhf4gZpbsIzG5\n9FI4kvkd8e+NnWy2AgMesZJ5F+A6zCT4wzG58TWLyHyEqj6YJP5ieJmXIE9C12/FaerUldDvufFE\nJ/W4E/r9wGm8dNoXdNsgYWVpCL2EwFXhym/gD6sSuxhGYwv9EcwEaUdRWuiOZD7Usmyf8jKPvNcZ\n+Ieq7hRdzCdIs/RX1VeSxF801lqLNM8yh9oKXb+N+ZnNgdDv+X1E4Fs6CAaSC/1x88/IyfDbYfDQ\nsbBrp4R1OhD65Lc2YLMNPi5bbtRceHwW9FwD9ivyx0TjCn0xplH7NUboyxUp50jm91uW/ZVVy/yn\nqvqFiCyFmR32KVW9M/QAtCvND0A3SUt4XuYWZCX0WHnHkbHQW8k7jloI/fHib41+D46+F+4/Grpv\nlDCmCoU++a34PwmKCf3uOXDvHPgF8D8L4U9Lw5vAkT+Bo2OG5ze20B8DvgB+A8QtFuJI5paJDjms\nVW+WoZj1GtbETNXZHzNJfmGdhgc0NPOLiFyE6Zq4EDhLVUv81CbDy9yStIRuLfAoKQvdSuBRshB6\nBb8KT02BI4fAvUfBnpskjKmM0IsJPEqc0FVNq/zZb+DSRdCvHXRfDfZrH59ugUYX+gjgP8DRQPvI\n+45kPtyybG8/aMg+gDqRObgTuu7V5KQe10KvSuBR0hB6grbMsx/C4XfD3UfCvp0TxhQRuq3Ao8QJ\nvZBikXaweBH0KpFqKdC4QldgFGZp4aOhxSIVjmQ+wrLs/l7m9gHUkcyheqEP1M9a7J+1t6N1OhII\nfVCkBXzS75KFsoSkQi/MPuhoRfsXPoJDB8OdR8D+CQaKTroARJM+VTVEhX7rbPjZMrBvexj9A0xb\nACetYldXY0pdgdHAVOBYmoXuSOaWDQTZz8vcPoA6kznYCT0q7zhqIfSowKPUTOjFpo91JPSXP4be\nd8Fth8NBW9ifNylmEYO0hF4J8tK85p0OlzmIJo8o8CRmDvtjMalpRzJ/xrJsdy9z+wDqUOYQL3Qb\ngUfJQujlBB4lM6Hbzv/tSOivTYMD74SbDoVDS8QWJ/AotRB6C4FHaWihPw1Mwgj9Kjcyf96y7G5e\n5vYB1JHMh7Nfi/2pCefTLuBM6LBE6pUKPEpqQq92AQdHQn/9U9j/drj+YPhVMPLfRt5xpC30kvKO\no2GFDvAMZmGSWW5kbjmgXrp5mdsHkHOZRwUeJS9CH/Bk8+sOCWMp4Ezosx3V40job34Ge19rZnM5\nIGFdroVescCjNLTQnwOedCPzcZZld/Aytw8ghzIvJ/AotRJ6WOBRvNBbc3WoBf4f4C7gXKB3wpCS\nCn2zhwKR75RQ5AUaWuiOcuaWUynINl7m9gHkROaVCjxKVkIvJfAoXugtBR5lBma43hmYmbWTUKnQ\nCwJvVY8XehkcyXyiZdmtvMztA6jVgs7vNb9+pHMykRdIQ+iVyDuOtib0UvKOYyamhf47zESsSSgl\n9GLyjq3HC70EjmQ+ybLs5l7m9gFkKPOwwKPkSejfiMMHojS+0CsVeJQvMS30EzGzgiQhLPRKBN6q\nHldChwaTuiOZf2hZdiMvc/sAUpZ5KYFHqaXQXQs8SrVCPymaUF4vaSQBCYU+KZj1bmTySACYhWmh\nHwsck6CelYEOjiY49UKPw5HMP7Es+7NsZB5Ml1uO2ar6p6J1NJrMK5F3HFkKPW2BR7EVeiuBR6mB\n0CeVmLLUldC/AgZjlkk4wfKclYscr7nQG0beURzJ3HJIiKyXmcw/Bi4GBNO5vlUR4EJVLTqGORcy\nNx9Q9XN+L3FiDychpSb0rOUdR5zQy8o7jgyEXkrgUVwJ/RtMC/0w4JQiZYoJPErmQm9YgYdxJPOZ\nlmXXykzmfVX1miRlciNzoCKhF/Vij8QhAe6EPkFSW4y7ajpQpcCjpCD0SgQexZXQZwP/wPRBPy04\nZivwKKkLvU0IPIwbmS+eZVd2qTV8ztw+gJDMobTQrRu2PRKFtIRqhd67f0uBD7jERTTJufD75tfL\n9XFUqQuh3w+TLH+5bHAh9W+BoUAv4DzM37jV4lTobU7eUdzIfIFlim+ZVbKVuYhshOkt24nmVa5V\nVcs2v3Inc2gWeqKsRI8E54awEXpU3nHUQuhhecdRM6EXWeWl1kKPpqBmA9cAewEXUjuhv3Do9kte\n7yYHJYiiEXAj83llfjcKLL9i5jJ/CxgEvI2Z3B2MzMtODZZLmQNws4OFkHskrwLihW4j8ChZCL2c\nwKNkJnTLZbqyFnq5h8LfAtdiplf/E9kJPSzwKG1b6G5kPnthuYWjDass/WPWMn9VVbtWdW4uZe5C\n5AV6uKtK70n+maYh9EoFHiU1oVsKPEraQq+0m+Z3mKXXd8T8pKYl9FICj9K2hD4LeDfYPnYi8680\nuoJRPKvLD1nL/BhgY8ySLPMLx1X19bLn5lLmkBuhH9P51iWvB/cv1r+hMpIKvX9oYv35jla1dyL0\nFTHrCTjApdDfclDH98D1mHU6/wIslaCugtArkXccjSv0xZhVht7FzGU+D+gCdAaGOpH5F7qSVdmf\nyndZy/wyzFCHKTSnWVDVPcqem1uZQ02EHpZ3HLUSev8SK6PUVOgrxhzLgdA7RNYCHulg/M1c4EZg\nM+BSqhN6YRrt9TWZyAs0jtDnAx9gBP4+8BOMwLsA69L8v+0mzfIftVvGaV2ZHV3Q+XZMR6eZqrpV\ncKwrcAOwDGbh5tNU9bXgvX6YYQuLgDNVteRvh4h8AGyuqj9WfF+5ljlkIvRyAo+SldBLCTxKpkKP\nE3iUGgg9KvAoLoT+A3ATpqvBFUA7i3OKrYPghT4b0/J+F/gEWJ9mga9a5Bw3Mp+ua1iVXV9mRWW+\nOybzNjgk87HApar6uIj0As5X1T1EZAtgCCZD1wEYA3RW1cXR64Tqfxg4VVVnVHxf5WQuIj2BgZif\n20GqenlMmeswvbjmAsep6oTg+KqYJ7M/x4xqOkEj08KXlTmkIvRKBR4lLaFXIvAoqQrdRuBRMhB6\nOYFHcSH0ecDNmEcEV9Na6JaL2ABtTeiKmXy4IPBvgE0x6ZNNAJsP043MP9a1rMpuIDNbXU9EOgGP\nhmQ+FHhIVYeJSB/gAFU9OmiVLy44U0RGAU1RB0bqfgbYGniN5px58q6JItIO87++N/BpcIE+Gppz\nTET2B05X1f1FZCfgWlXtFrx3F/CMqt4uIksDK6rq7Mg1ysu8QAKp9zql+cnTmrhJyLoQet8Bl3LN\n6H4OojE4E/pJbupxLfRK5R2HC6HPx7RSVsf0drFcuCaWxhb6QsyizIX899I0t747Yve3TRg3Mv9Q\n17Uqu5H8x0bmG2C+wxWTD9pZVacF8628rKr3BOUGASNV9YESsfWIOZy8a6KI7Az0V9Wewf6FQc2X\nhcrcBDzWSvMxAAAgAElEQVStqvcF+5OB7pgGzARV3ahkAJXIHCoSeljgUWop9L4DLm11LA9CXy7a\npLwzaSQBCYU+MviodnUg8iV1OhD6j8BVQHvgdJpHeFSDK6FDHqT+PSbv/S7wIbAWzQJfk2T9gdzI\nfFKR6YpfHTuPV8c2/3DcOGC2jczHADeq6kMi8ivgFFXdp4jMR6iqo2FkLSkn88OB/VT15GD/aGAn\nVT0jVOZRTL7oxdCNXYBJ+N8MvIPpBDAeOEtV50auUZnMoaTQSwk8SpZCjxN4lFoIvZXAo9yZNJKA\nCoU+sshHk0ehX4MR+VnUTui1F/iXNHcfnAFshEmfdKa6HF0x3Mj8bd3YquyW8oGNzL9V1ZWD1wJ8\no6qrRBu/QZqlv6q+EhPTY6p6YJm4S5YpJ/NfAj0tZH6Zqr4Q7I8Bzsf8ufESsIuqviYiA4FvVfXi\nyDUqlzm0EHolAo+SltBt5B1H2kIvK+847kwaSUAJoReTdxx5E/pCzEMlBc7GdGmoFluh117ei2ju\nPvgusIDm1ncnkv0vlMKNzN/UzlZlfyHv2cj8daCvqj4jInthfLhj6AFoV5ofgG4SN02siMwGni0T\nzpaqumHR+yoj826YhH0hzdIioR8cuwkYq6r3BvuFNIsALxUuLiK7YaZwPDByDTXFC3QCisYLwPa6\nNwBrUfED31hcCX0NR/W4FnpVAo9yp4M6oIXQKxF4lDwK/TqM0voCduML4ykm9NoLfB4tuw+uQsvu\ng2l0x54KfBTaf8aJzF/Xza3KbieTor1ZhmKEtSbmT5CLgYmYXqvLYTo8nRbqBHIRpmviQkxmIraL\nQ5FceZT5qvpSsTfLyXxpzCe3F/AZ8CqlH4B2AwaGHoA+C5ykqu+JSBPQXlUviFzDqmVeEHiUWgt9\nHDu02N+PBN1RQiQV+hP7NjfL9xntwua4Efr3MHK4g3rIp9BvxGSMz8ON0Gsv8G9o7n0yDfPQsjCA\np1j3wTRx0zJ/Vbe0KttV3m6cWRODfpOFrom3qeqlInIqgKreHJS5AeiJ+Tk+vjD0VER+gXnovyzm\nK/34SnqzFBN4lKyFHhV4lFoJPSzwKDUTepGpBlwJHaqT+sqHtj52X4LpdwssAv4OfA38AbvOdmF+\nHVo1XrZpSh5QxSymZffB2Rhxd8GMMl+uBjGFcSPzl9RuVfCd5Y3GkXnqAYRkbivvONIUejl5x5GF\n0EvJO47MhG45V0zWQo8TeBQXQl+MefI/E9MToFxoYYFHyUboC2jZfXBZmtMn61N598E0cSPz562f\nTYz3MrcOQES3VzeScSn0agQeJQ2hVyrwKKkJvcrJvtIWuo3Ao7gS+iDM4IwLgBUi75cSeJR0hP4d\nzd0HpwJr07L7YF5xI/OndGersnvKS1nPzdIbeKzUKNGi5zaSzCGZ0EfKxCWvNy8/4MoKF0IfwlHc\nzW8cRGNwJvQb3VTjUugAvRxMGuZK6HdgHt/1A46vQOBRkgtdgS9oTp/MxKRNOmNGYbrsPpgmbmQ+\nWu0aRfvK81nL/B5gZ+CfwO2qOtn63DzIvP3sr9hi5Xec1Wkr9LC846iV0IdwVNH3XEm9KqG/H3PM\n0QjPJELvFf2YHHkpqdB//Tiowll3wUtvwuibYLVq15+jGqEvwsx5Uug+uIiW3QeT9IqvFW5kPlJ7\nWJXtJWMzT7OIyCpAH+A4zLfwHcBQVZ1T8ry8yBzIROjlBB4lK6GXEniUTIUeJ/AoNRB6K4FHqZHQ\nfx3zMavCuXfD2HHwxE2wRoJOIOWFPg8ze+q7wb+r0izwdUin+2CWuJH5Y7qXVdkD5cma5MxFZE3M\nVLhnYwZebgpcp6rXFT0nTzKHdIReqcCjpCX0SgQeJVWh2wg8SgZCLyvwKBkJPU7gUVThwnth1Isw\n5ib46erVx9Na6F/TnD6ZDmxA8+hLu6le6wc3Mn9E97Uqe7CMzjrNcjCmRb4pMBi4U1VnisgKwDuq\n2qnouXmTObgR+ngpjJgt+ZeJNa6EPgu7qTdtcCb0Gx3l0B0LvWJ5x5GC0G3kHYcq/Pmf8PDT8OQt\nsHaVPwqLF0O77U6iuffJHJq7D25E7bsPpokbmT+gvazK/lJGZi3zuzBdwFuNBhWRvVV1TNFz8yhz\nqFzozfKOo3ZCnyTjWx1by7JVYEO1Qv8h0r+i942OTJy0mkL33wQPC1vhQujb0BxbAlThkofh3sfh\nqVtg3Z/anffDPHjyVRg+Fh591uTee3eHy+84AdN9MMnaR/WEG5kP1UOsyvaRhzOTeTBIc4yqZUI/\nen5eZV6glNRLCzxKdkKPE3iUWgg9KvAoNRN6MUnWWuhxcTkQOsBfh8Pgx4zQO6wdX2bGLHjsWRj+\nDDz9Gmy3GfTuAededQY4/AuvvnAj87v1l1Zlj5YHsm6ZPwn8UlW/qfjcvMscWgq9MoFHSU/oNgKP\nkoXQywk8SmZCt5VilkK3jcmR0K8YAbc8aIT+s3VNq/3fHzS3vid9CPvtAsNGH4ZZvKGyz7IxcSPz\nO/UIq7LHybCsZT4c2BbzG1SYYVZV9cyy5+Zd5j+sUnh4+xNHV3Qj9MrXeI/HtdArlXccqQg9iQDT\nFHq1cTkS+jWj4W+DYb+dYex4kw/v3QOuH3oM5kFmPXYfTBM3Mr9N7TofnChDspb5bzFdjgpiFozM\n7yp7bh5l3izwKF7ocRzFkCWv9+DppOEAjoQ+FXdjUVwKfVdH9TgS+jHXwlLtYPCjv8OMxKz37oNp\n4kbmt+gxVmVPkX9kLfOzVXVguWOx5+ZF5sUFHsULHVoKPEpNhT415lgehD4osn9HkkBCJBS6fBn8\n/v1mQPJY2gRuZP53Pc6q7O/lzqxlPkFVt40ce0O1/MxguZB55YtTtD2hl5J3HJkJPU7ecWQt9Ki8\n46iB0JfIOw4vdAvcyPwGPdGq7OlyWyYyDxaCPgrYHXgu9NZPgEWq5Uc51anMoS0IvVKBR0lN6LYC\nj+Jy+o84qdsIPEoGQi8p8Che6GVwI/Or9DSrsufK/2Ul8w0wq/JchpmbrXDNb4G3VHVh2TrqV+aQ\nTOiuJB7GhdA/5Wx1N+jDmdDPc/RQ1LXQqxF4lBSEXpHAo3ihl8CNzK9oXvmyJOfL9bWYm6UTZnm5\nMcHIz3bl5mWBupc5VCb0NAQepRqhf9rqSB6EfiXntdh/7jxHPW+SCP10NyG0IqHQvz6/PQCrD5lb\npqQlXuhFcCPzS/Vsq7L9ZGDWOfNTgJOB1VV1YxHpDPy9wdMsYYoJPQt5x1FO6K3lXQxXUrcRelTe\ncdRE6GkJPEqFQi8IPIoXepq4kflf9Fyrsn+Wq7KW+ZuYBaBfLjwIFZGJhcWjS9EgY4DnRF7PoXYi\nh3hZfxra7Bko810ExNPsEXv8Ss5bstmw+98cpVvKLWZxemjLiuPLF/n6/PZLtmJ8dZSjwT339HdT\nj6cVi1jaaosiIreLyAyR5tn7ROReEZkQbFNFZELovX4i8r6ITBYRm5bQfFVd8ksfDPG3anE3SMu8\n8XGZdhmH3ZJZpUilhZ6luEsRaaGXEncpfAs9Ddy0zP+sF1mV/Yv8b4vricjumGWaBse1lkXkb8A3\nqvo/IrIFMATYEfPn+higc6lVhETkSswq2sdifiNOw8yW+MdysfrhZXXCQJlftdAP558t9l3IfPe/\njU4u9AsSh5EO38PXA6oTeJivjlrBjdALLXQvdWcsqnJdU1V9LnhA2QoREeAIWPJn8MGYRSUWAB+J\nyBSCFEqJS1wInAhMBE4FRmD5mN/LvI6oROhRgYe5jH5cyKWJ46lK6HkV+A0td1fr/0M+hO4Fngrz\nWTaNancHZqjqB8H+erQU93TKPFBT1UXALcFWEV7mdUYpoZcSeJRMhV4nAo9SM6F7gadOXD4cYNrY\nD5k+9sNqq+0DZQeHlMxri8hBwCW0XNdPVbXsgoM+Z16nFIReicDjcCF0iOTQ61TexXAhdCiRQ/fy\nrgA3OfMz9AqrstfL+a2uF6RZHg3nzIMHldOB7VT1s+DYhQCqelmwPwror6pFp34VkQ+AQ4G3S+XW\n4/At83rjKpM/HXg1PH/ODomrc9VCn38JLFeuh0otqFLgYVJpoXuB15Rqc+Yl2BuYVBB5wHBgiIhc\njUmvbAq8Wqae6cC/KxU5eJnXB1fFd1Hb7epxNRX6mLkt0yvzV8yJ0MPLum1KdWuaRnAh9NWfnWsW\nBeruRV5rFlYpcxEZCnQH1hCRacDFqnoH8GugxWqxqvqOiAzDLMi8EDhNy6dCLgBGisjTwI/NVenV\nZWPzaZacUkTgcbgQOtilXKICj6MmQi+3LqcDoUNlKZfVny2RJ/dCrxI3aZYT1e5Pttvk9KwHDT2B\nGSQzEVjSOlfVsj8wXuZ5oQJ5x5GW0G3kHUfqQq9mUeUMhF5S4FG80KvAjcyP079blb1Tfp+1zN9W\n1S2rOtfLvJb0havKPqS2xpXQd5hb+RJ4cTgXejUCj5KC0CsSeBQv9ApxI/Oj9DarskPkxKxlfgXw\npKpW/NPeIMP564m+oQ0491tnNe929biqz91h7vglmyvmu5gh8b3Q5oJN3VRzDX1Z/dm5yUQO8Iwf\nsl8LFtLOaqsBp2Fy5vNEZE6wWUnCt8xTp69dsYxb6C6lXY6KWujlpF3tXOpRKmyhXzwgfvj3X579\nq4Ng8C10a9y0zH+pd1uVfUCOznwK3GrxMk8FS4FHSVnoWQo8SkmhV9rqzkjoxQQexQs9S9zI/BAd\nWr4g8LD0yTrN8iBwGzCy0u6JdSjzwp+lefvBr1LgcTiS+vPn7FBTgUdpIfSkaZOUhG4r8Che6Fnh\nRuYH6TCrso/KEVnLfB/M/J3dgGHAHar6rtW59SHzYnnFWv/gbxbZ7+Wu6iRCDz2on/dmfv5CfHyF\nfen9gqMpdMGJ0OWYJv6sP5YvaIEzoYOXelHcyLyXPmBVdqT8siZpFhFZFTgS+BPwCXArcHcwaVf8\nOfmUeSUPhbL8oY/KO44aCL1ML6taCf3xFeK7NdZa6HJMU6tjuRC6F7gFbmS+rz5iVXa0HFyLZePW\nAI4BjgY+w8z3shuwpar2KHpefmSe5Kl+mr8ENgKPkoHQ7brJLiEroRcTeJSshR4n8Cg1EboXeIW4\nkXkPHWlVdqz0yjrN8hBGOv/ApFj+E3pvvKoWnb86JzJ3EYPLX4pqBB4lBaFXKPAoaQndVuBR0ha6\njcCjZCJ0L/AEuJH5bmr3s/e87Ju1zPdX1RGRY8uFVx8qem45mYtIT2Ag0A4YpKqXx5S5DmOvucBx\nqhpeNqkdMA6YrqoHxZzrSObgTuguZA7JhX5N88tN3PRHdiH05Vd8k0f0Dw6icSt02a3JST3Ohe7l\n7RA3Mt9Zn7Iq+5LsmbXMJxTW/gwde11Vtyt7bimZByJ+FzMj2KfAa0AfVZ0UKrM/cLqq7i8iOwHX\nqmq30PvnANsDP1HV3jHXUFz6vO6Ffk3xt2oo9OVXfLPVsTwI3ZXAoyQV+l8klcUPPI5k3lWfsSr7\nqnTPROYisi5mMYt7gKMAwVhxZeAmVS0rpHIy3xkz/27PYL/F/LzBsZuAp1X1vmB/MtBdVWeIyPrA\nncBfgXOKtswliKHNCr2EwKNkKPQ4gUephdDTEniUSoXuBZ4FbmS+vT5vVXa87JaVzH8LHAfsgMlk\nFJgD3KmqD5aro9wUuB2AaaH96cBOFmU6ADMwlvoD5tulPIXvosT0x43QJ+NG6CNpLfQKBB5mygAn\nQl/+F9pK6DbyjnKwXOlE6MN3Ld5tMSt5R/mLLFtS6F7e9UsK85knQlXvAu4SkcNVtaoVZ8rJ3Fat\n0W8uEZEDgZmqOkFEelhHlFToS7JN/WFCnoRepbzjcCh0pryVuJ40hF4rgUeJCt0LvDGYT3WLo2fA\nYyLyG8yyce0IjKiql5Q7sZzMPwU6hvY7YlrepcqsHxz7JdA7yKkvD6wsIoNV9dhWV9Gm0E4PkB6V\nCX3bYsdrLfTJDq5dhCRCnxL9Pzk0cTguhH6w7ALskjgW13iB15KpwEfOa81byzzEI8A3wHhgXiUn\nlsuZL415ALoXpvP6q5R+ANoNGBh+ABqU6Q6cVzZnHqWU0IsJPA4nQgc7oaco8Dhshd5K4FGSCx0q\nz6EbgXs8trjJmW+sb1uV/UC2rJv5zEu2zFV1oYicjplJuh1wm6pOEpFTg/dvVtURIrK/iEwBvsfM\nKxBbXcXRhVMulcg7Sqot9IzlHaVYC72svKM8RBYtdC9vTx6o0fS2NrwoIlurasX5z3wMGirWMj8x\n+NfVXFHOWug5ZJP+VQg8DvctdC9wjzvctMzXV7s5kKfLpi2uJyK3AwdgngduFTp+BmYu8kXAv1T1\nguB4P+CE4PiZqvGjlURkYvCyHWbW/alAYaCQqurWZe8rdzI/sUjBmgs9OiuiwweaucOF0B9yUIfH\nE8WNzNfVD63K/kc2isp8d+A7YHBB5iKyB3ARsL+qLhCRn6rqFyKyBWZelR0xPfzGAJ3jprYVkU6l\n4lDVj8rFWu4BaDYUE3iY7XEj9IpSLqWmte1L4wq9mpSLl7enfqj2AaiqPhcj3t8DlxZmNFTVL4Lj\nBwNDg+MfBanorsDLMfV+FN4XkbUwHUesyYfMbUld6JXOSd7Whe4F7qlP5v/otIfSpsB/icj/Ynqg\nnKeq4zAjOsPiLozBKYqI9AauCs6dCWwATAJ+Xi6I+pI5uBU6wISka3C2NaF7gXvqn0ULnapvaWA1\nVe0mIjtiFpXYqEjZcnnt/wF2Bp5Q1W2DFM4xtkHUH0mFfkjo9YSipSqg0YXu8TQWixbGp1kWPfc8\ni5+3G+ofYjrwIICqviYii0VkTYqPwSnFAlX9UkSWEpF2qvq0iFxrE0Q+HoCeVGUMtkI/pMz7/ZO2\nzgs0qtA9nrzg5gHosrNmW5X9cY1VWl0vyJk/GnoAeiqwnqr2F5HOwBhV/VnoAWhXmh+AbqIlpCsi\nYzB/Dl8KrIlJteygqmW7hC1ldUd5peg07RiBF7ZyDHC1kLLDdUA9Hk9qLFzQzmqLIiJDgReBziIy\nTUSOB24HNgq6Fw4FjgVQ1XcwKZd3MBM0nVZK5AGHYKYS7wuMAqYArQZbxlHfLfMChRa6jbhL4Vvo\nHk/OcdMy51PLkfIdlq/JGqAAIvJTYFZcV8Y46rtlDmbM6WYkFzm4a6EPcDNNrcfjSYmF7ey2jBCR\nnUVkrIg8KCLbicjbwETgcxGxWhSh/lrm35d4z9X049W00OO+CPo38IhTj6cmOGqZv2vpnC6S1Xzm\n44F+wCrArUBPVX1ZRDYD7lXVbcrWkQuZb6ul89+lBB4lS6HbtOS90D0ehziS+b8tvffzzGT+RkHY\nIjJJVTcPvddqKbk48pNmifZM+T60VYKrea+KiXrAys2bVT0+5eLx5I6Fllt2hL9dKpr6tkB+WuYF\nXLWsXdXjCt9C93gc4KhlPt7Se9tn1jJfhOnFAtAe+CH0dntVLTsmKH8yh9oLfVpkv2NsqcpJLPRw\n10ffY8bTFnEk85ctvdctG5m7ID9pljCuUiWV1DMttMW954KqUi59Q1v0uMfjqYpFllsdkc+WeYE0\nW+jVCNpVCx1KtNIrlbRvoXvaEo5a5s9Yeq97/bTM8y1zcCt0Fy3sVISetJXthe5pKziS+eOW3tuv\nfmSezzRLmKQpl8nB9rCDWMBdymUy0Kc/btIlPuXi8VRE/nqzJCb/MofKhT45tIVxtVpRtUKPi6uP\nnxfG48kcL/MaUkrokyku8ChZCt02Li90jydbvMxrzOTIaxt5x5Gm0KuNywvd48mOBZZbHVF/i1O4\n6rY4ntJTCNgyjcpHqRajz8ow1MXMjQWh+wejHk8sddbt0Ib6k7lLqhX62JhjOyYLZQmJhP6CoyA8\nnganzlIoNrRtmYO90MeWef81aiT0UgLvCrzqICCPp8GoavaTfONlDsWFPrbCejITeiUtcC90j6cV\nvmXewBQeis5JWE8qQk+aPvFC93ha4GXegEwY2XJ/E6tFPUrjQuhDw8P9uyasrFCHF7rHA3iZNwxR\ngYeZMrJ2Qh9abL6WV/FC93gcUmW3QxG5HTgAmKmqWwXHmoCTgC+CYhep6sjgvX7ACZj+M2eq6uhE\ncZcgHzKfAmyS8jVKCTxKlkIvKvAoXugejzOq75p4B3A9MDh0TIGrVfXqcEER2QL4NbAF0AEYIyKd\nbRdorpR8yBzcC70SeceRltCt5R2HF7rH44Qqe7Oo6nMi0inmrbjJuA4GhqrqAuAjEZmC+eV7ubqr\nlyY/MofkQp/gYsBNCFdCTyTwKF7oHk9i3OfMzxCRY4FxwLmq+g2wHi3FPR3TQk+FfMkcKhd6UYHv\nipNBNNUKfUqay8QlEboXuMfjeKj+34FLgtd/Aa4CTixSNrU5x/MncygvdOsWeIZCT1XecdgK3cvb\nkx6qzatnidTROrfFcubTxsL0sRVVpaozC69FZBDwaLD7KS1XQFg/OJYK+Vic4idFYigIPXH6xNEw\n96jQMxd4HHFC9wL3pEdY4FHSF7qjxSnOsPTe9a0Xpwhy5o+GerOsq6r/CV73BXZU1aOCB6BDML+k\nHYAxwCaaknTz2TIHmDMIJgAc4aAyhy303Ikyb/F4GpFSAo+Wq4sWepU5cxEZCnQH1hSRaUB/oIeI\nbINJoUwFTgVQ1XdEZBjwTnDF09ISOeStZT5nUJFSLoQO1QvdC9PTtrCVdzHSE7qjlvlJlt4bVD/L\nxlm1zEWkJzAQaAcMUtXLY8pcB/QC5gLHqeoEEemI6Y+5FuZb6xZVva7VBYpKvMAwsm2he3l72h5J\nBR6tK9ct9Pm1DsA9ZWUuIu2AG4C9Mcn710RkuKpOCpXZH5ML2lREdsI83e2GeWbcV1XfEJGVgPEi\n8kT4XHvSFroXuKft4VLgcXXnVuhtdDh/V2CKqn4EICL3YjrDh4XcG7gLQFVfEZFVRWRtVf0c+Dw4\n/p2ITML0vaxC5uBW6OAXb/C0NdKUd7Hr5VLodbaKkA02Mu9AywXSpgM7WZRZH5hROBA8Ad4WeKWK\nOEMkEbpfvMHT9sha4HHXz53Q2+hKQ7ZPSKMPCZacF6RY/gmcparftT51eOh1l2ArRSVC94s3eNKj\nIMq8yarWAo9SvdCnAh85joY2m2aJdnzviGl5lyqzpHO8iCwDPADcraoPx1+it120LSgldL94gyc9\n4kSZh9Zn3gQepbr/ow2DrcAzboJpozIfB2wapEk+w8wC1idSZjhwOnCviHQDvlHVGSIiwG3AO6o6\n0FnUSwgLPUkKxQvdU5y89rHOu7zjyMOXHtCQOXOrfuYi0ovmrom3qeqlIlLoGH9zUOYGoCdmrfrj\nVfV1EdkNeBZ4i+a0Sz9VHRWqW+HWKkIPjwrdvIrz4/BC9xiSiDJNWdWjwOOo/v/IUT/z3S2zx8/V\nTz/zfAwaspZ5qWH9XuieZLgUpUuhN4rAo1T3f+RI5jtbeu8lL3P7AMrKvJJ5WbzQPfakLclqhd6o\n8i5GZf9PjmS+g6X3xnmZ2wcgotBE8zS/SSfV8kL3FCdrUdqKqq0JPIq90B3JfFtL703wMrcPYInM\nAX7iqFYvdE8ztRZlMVHVOq68YSd0RzLfytJ7E+tH5jmbNXEOboQ+ieqF3inyeljSYDw1JA/CDPfg\nyEM8Hhqya2LOWuYFsm6hdyrzvhd6PeMFml8qf67gqGW+iaX3pviWeULSbqF3qrCeI/BCr19EBnih\n54hc9DNvi7Mm1g6XQgczO28SvNDrGS/02pELeUdpwDRLjmUOyYU+x1UgAV7o9YwXenbkUuBhGnAE\n6FK1DqA8lQh5TmQL40rCrlY98tSC3EvGkw2LLLc6og5kDqWFXkzecXihe7zQsyD3fwEttNzqiDqR\nObSUdSUCj+KF7vFCz4JcC71KmYvI7SIyQ0Qmxrx3rogsFpHVQ8f6icj7IjJZRPZN5V4C6kjmUL3A\no/i8t8cLPQtyK/QFlltr7sBMKNiCYL3jfYCPQ8e2wMwyu0Vwzv+JSGrOrTOZu6RaoQ8LbZ56RbV/\nfkXTYOTy/7nKlrmqPgd8HVPj1cD5kWMHA0NVdUGw7OYUzHzbqZDz3ixpY7tikRd3I5BLqbQRcjOP\neQqIyMHAdFV9yyzhsIT1gJdD+9NpnoTKOW1c5lBc6F7g9Y6Xd77I6xJ7SRCRFYCLMCmWJYdLnJLa\nkHsvc8CLu3HwAs8/+W6ljw02azbGDCl/M2iVrw+MF5GdKLGcZhrkdG4Wj8ceL/D6pKaLU/CjZell\nW10vWELzUVXdKqbuqcD2qvpV8AB0CCZP3gEYA2yiKUnXt8w9dYkXeP1T2xZ6dZ3IRWQo0B1YQ0Sm\nARer6h2hIktErarviMgw4J3ggqelJXLwLXNPneDl3bjUZKUhZluWXsXPmujxJMULvG1Qmxb6Dxlf\nL328zD25wgu8bZK90Btvpi0vc09N8fL2FMhW6HU28YoFXuaezPEC9xQjO6H7lrnHUxVe4B5bshG6\nb5l7PFZ4eXuSkL7Qfcvc4ymKF7infvC9WTyeFniBe1yRbW8Wn2bxeLzAPc6o3QhQn2bxtEG8vD0u\nycckW75l7mkjeIF7XJIPgYfxLXNPA+MF7nFF/uQdxbfMPQ2El7fHJfkXeBjfMvfUOV7gHpfUl8DD\n+K6Jnjqm1iKP+8WvdUyeyqhfeUfxLXNPHSMyIFN52vziZx2Tx2NovJz5UuUKiEhPEZksIu+LyAVF\nylwXvP+miGxbybmebEm7ZSUyYMlWyTme+qBxvngXWG71Q0mZi0g74AagJ7AF0EdENo+U2R+zrt2m\nwCnA323PbWym1jqAolQrz7Fjxxatr1KBu4rJBcXuq95J674aQ+gLLbf6oVzLvCswRVU/UtUFwL3A\nwZEyvYG7AFT1FWBVEVnH8twG5qNaB1CSauRZkENY3i4lXCuhe5lXTv0LvY21zDErSk8L7U8PjtmU\nWc/iXE8NqTQVMmDA2EzSNJ76oL6F3vZa5rarPdfFgqee1pSSZxqtbxu80OuH+hX6D5Zb/SCqxX0t\nIpz1k4UAAAaOSURBVN2AJlXtGez3Axar6uWhMjcBY1X13mB/MtAd2LDcucFx2y8Mj8fjQVUTNR4r\ndU7S62VFua6J44BNRaQT8Bnwa6BPpMxw4HTg3kD+36jqDBGZZXFu3fxHeTyexqBRnVNS5qq6UERO\nBx4H2gG3qeokETk1eP9mVR0hIvuLyBTge+D4UuemeTMej8fTVimZZvF4PB5PfVB20JALRGR1EXlC\nRN4TkdEismqRcreLyAwRmRg5fqWITAoGJT0oIqtkEXcpHNyT1flZU8F9xQ4IE5GuIvKqiEwQkddE\nZMfsoi9O0vsK3jsj+Dl8W0Qujzs/S1zcU/D+uSKyWERWTz/q8jj4GcydL7IgE5kDFwJPqGpn4Mlg\nP447MIOMoowGfq6qvwDeA/qlEmVlJL0n2/OzpmxcZQaEXQH8WVW3BS4O9vNAovsSkT0wYyq2VtUt\ngb9lFXgJkn5WiEhHYB/g40witiPpfeXRF+mjqqlvwGRg7eD1OsDkEmU7ARNLvH8ocHcWcad5T5Wc\nn7f7AnYGRoX2LwQuDF4PBY4IXvfJw2fl6L6GAXvW+j5c3lOwfz+wNWbI8uq1vidX9xU6ngtfZLFl\n1TJfW1VnBK9nAGsnqOsEYETykBKT9J5c/p+4xCauUoPJLgSuEpFPgCvJT6so6X1tCvyXiLwsImNF\nZIf0QrUm0T2JyMHAdFV9K9UoKyfpZxUmL75IHWezJorIE5hv0Sh/DO+oqlbbt1xE/gj8qKpDqjm/\niuulfk8uzq8UB/dVKtbbgDNV9SER+RVwO+bP+NRJ+b6WBlZT1W7Bc4BhwEZVB2tJWvckIu2Bi2j5\n2WTWZS/lz6pwjUx9UWucyVxVi/7CBg8A11HVz0VkXWBmpfWLyHHA/sBe1UdZGSnfU+L/k2pxcF+f\nAh1D+x0xLSOArqq6d/D6n8AgFzHbkPJ9TQceDK7zWvDAcA1VneUo/FhSvKeNMem/N0UEYH1gvIh0\nVdXUfxZT/qxq4otak1WaZTjw2+D1b4GHKzlZRHoCfwAOVtV5jmOrlkT35OD8tLCJa8lgMhFZFjMg\nbHjw3hQR6R683hPzACoPJL2vhzH3g4h0BpZNW+QWVH1Pqvq2qq6tqhuq6oYYEW6XhcgtSPRZ5dQX\n6ZNFYh5YHRiD+cUeDawaHF8P+Feo3FDMaNH5mHzY8cHx9zFP2ycE2//V+mGDg3uKPb/WWwX31Qt4\nF5gC9Asd3wF4BXgDeAnYttb35Oi+lgH+AUwExgM96v2eInV9SH4egCb9rHLniyw2P2jI4/F4GoCs\n0iwej8fjSREvc4/H42kAvMw9Ho+nAfAy93g8ngbAy9zjaWOIyK9E5N8iskhEtitSpqOIPB2Ue1tE\nzgy9FzuRlYgsKyJ3iMhbIvJGqIsqInK8iEwMzhkpImuUifE3Qdm3ROQFEdna1f03Kl7mHk/bYyJm\nzpJnS5RZAPRV1Z8D3YD/tpjI6mTMamJbY0aWXgVG8piJyboH57yFWdCmFB8C/xXU9Rfglspuse3h\nZe7xtDFUdbKqlhzMpaqfq+obwevvgEmYft6o6hOqujgo+gpm9CjA5sDTQZkvgG+COWwWAl8DK4kZ\nbroKZgQnIvJTEfmnmGmTXxWRXYLzX1LV2THX8BTBy9zj8ZREzNKP22KkGiU8kdWbQG8RaSciGwLb\nAx0D8Z8FvI2R+OaYOXwArgWuUdWuwOHET/9wIm1ksqwkOJubxePx5IcSE1ldpKqPVlDPSpg5ds4K\nWujh96ITWd2OEfU4zAjMF4FFIrIycB3wC1WdKiLXYyb5+iuwN7B5MD8MwE9EZAVVnRtcYw/MF8au\ntjG3VbzMPZ4GREtMZGWLiCwDPICZD/zhyHvHEZnISlUXAeeEyryAyalvDkxV1anBW/cDhZWBBNhJ\nVX+Muf7WwK1AT1X9Oun9NDo+zeLxtG1ip70Nctu3Ae+o6sDIe7ETWYlIexFZMXi9D7BAVSdjHmZu\nJiJrBkX3Ad4JXo8Gwj1ltgn+/RlmlsqjVXVK4rtsA/i5WTyeNoaIHIpJe6wJzAYmqGovEVkPuFVV\nDxCR3TC9Xd6iee7wfqo6SkTeB5YFvgqOv6SqpwW59VHAYswsjCeq6rTgmsdivgAWAx8Bx6nq10EX\nxRsxrfelgWeCugZhetx8ElxjQZBX9xTBy9zj8XgaAJ9m8Xg8ngbAy9zj8XgaAC9zj8fjaQC8zD0e\nj6cB8DL3eDyeBsDL3OPxeBoAL3OPx+NpALzMPR6PpwH4f8tiPr1bA1nWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f2bfbf8f990>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.pcolormesh(lons, lats,bathy,vmin=131,vmax=211)\n",
    "cbar = plt.colorbar()\n",
    "cbar.set_label('Bahtymetry [m]')\n",
    "plt.plot([lons[j-stp,i-stp], lons[j-stp,i+stp]], [lats[j-stp,i-stp], lats[j-stp,i+stp] ] ,'k')\n",
    "plt.plot([lons[j-stp,i-stp], lons[j+stp,i-stp]], [lats[j-stp,i-stp], lats[j+stp,i-stp] ] ,'k')\n",
    "plt.plot([lons[j+stp,i+stp], lons[j-stp,i+stp]], [lats[j+stp,i+stp], lats[j-stp,i+stp] ],'k')\n",
    "plt.plot([lons[j+stp,i+stp], lons[j+stp,i-stp]], [lats[j+stp,i+stp], lats[j+stp,i-stp] ],'k')\n",
    "plt.axis([-123.1,-123, 48.7,48.8])\n",
    "plt.plot(lon.values[0], lat.values[0],'r*')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The point at find_closest_model_point is too shallow but there are neighouring points at an appropriate depth. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}