working/interpolation_testing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test interpolation routines with some sample data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import numpy as np\n",
    "import marc_analysis as ma\n",
    "\n",
    "from ipywidgets import interact"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Open a sample dataset, and select a portion of its grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/daniel/anaconda/lib/python3.4/site-packages/xarray/conventions.py:382: RuntimeWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using dummy netCDF4.datetime objects instead, reason: dates out of range\n",
      "  result = decode_cf_datetime(example_value, units, calendar)\n",
      "/Users/daniel/anaconda/lib/python3.4/site-packages/xarray/conventions.py:401: RuntimeWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using dummy netCDF4.datetime objects instead, reason: dates out of range\n",
      "  calendar=self.calendar)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<xarray.Dataset>\n",
      "Dimensions:              (cosp_prs: 7, cosp_tau_modis: 6, ilev: 31, lat: 24, lev: 30, lon: 144, nbnd: 2, slat: 95, slon: 144)\n",
      "Coordinates:\n",
      "  * lev                  (lev) float64 3.643 7.595 14.36 24.61 38.27 54.6 ...\n",
      "  * ilev                 (ilev) float64 2.255 5.032 10.16 18.56 30.67 45.87 ...\n",
      "  * cosp_prs             (cosp_prs) float64 900.0 740.0 620.0 500.0 375.0 ...\n",
      "  * cosp_tau_modis       (cosp_tau_modis) float64 0.8 2.45 6.5 16.2 41.5 ...\n",
      "    time                 object    2-06-15 00:00:00\n",
      "  * lat                  (lat) float64 -21.79 -19.89 -18.0 -16.11 -14.21 ...\n",
      "  * lon                  (lon) float64 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 ...\n",
      "  * slat                 (slat) float64 -89.05 -87.16 -85.26 -83.37 -81.47 ...\n",
      "  * slon                 (slon) float64 -1.25 1.25 3.75 6.25 8.75 11.25 ...\n",
      "  * nbnd                 (nbnd) int64 0 1\n",
      "Data variables:\n",
      "    hyam                 (lev) float64 0.003643 0.007595 0.01436 0.02461 ...\n",
      "    hybm                 (lev) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...\n",
      "    hyai                 (ilev) float64 0.002255 0.005032 0.01016 0.01856 ...\n",
      "    hybi                 (ilev) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...\n",
      "    cosp_prs_bnds        (cosp_prs, nbnd) float64 1e+03 800.0 800.0 680.0 ...\n",
      "    cosp_tau_modis_bnds  (cosp_tau_modis, nbnd) float64 0.3 1.3 1.3 3.6 3.6 ...\n",
      "    P0                   float64 1e+05\n",
      "    date                 int32 20615\n",
      "    datesec              int32 0\n",
      "    w_stag               (slat) float64 0.0005467 0.00164 0.002731 0.003819 ...\n",
      "    time_bnds            (nbnd) float64 529.9 530.0\n",
      "    date_written         |S8 b'11/03/15'\n",
      "    time_written         |S8 b'16:54:31'\n",
      "    ntrm                 int32 1\n",
      "    ntrn                 int32 1\n",
      "    ntrk                 int32 1\n",
      "    ndbase               int32 0\n",
      "    nsbase               int32 0\n",
      "    nbdate               int32 10101\n",
      "    nbsec                int32 0\n",
      "    mdt                  int32 1800\n",
      "    nlon                 (lat) int32 144 144 144 144 144 144 144 144 144 144 ...\n",
      "    wnummax              (lat) int32 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...\n",
      "    gw                   (lat) float64 0.03071 0.03109 0.03145 0.03177 ...\n",
      "    ndcur                int32 530\n",
      "    nscur                int32 0\n",
      "    co2vmr               float64 0.000367\n",
      "    ch4vmr               float64 1.76e-06\n",
      "    n2ovmr               float64 3.16e-07\n",
      "    f11vmr               float64 6.534e-10\n",
      "    f12vmr               float64 5.35e-10\n",
      "    sol_tsi              float64 -1.0\n",
      "    nsteph               int32 25440\n",
      "    AREL                 (lev, lat, lon) float32 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...\n",
      "    AWNC                 (lev, lat, lon) float32 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...\n",
      "    CLOUD                (lev, lat, lon) float32 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...\n",
      "    EFFLIQ               (lev, lat, lon) float32 10.0 10.0 10.0 10.0 10.0 ...\n",
      "    FREQL                (lev, lat, lon) float32 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...\n",
      "    ICWNC                (lev, lat, lon) float32 5.41217e-11 5.41331e-11 ...\n",
      "    LIQ_ICLD_VISTAU      (lev, lat, lon) float64 nan nan nan nan nan nan nan ...\n",
      "    LWC                  (lev, lat, lon) float32 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...\n",
      "    PRCO                 (lev, lat, lon) float32 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...\n",
      "    Q                    (lev, lat, lon) float32 1.813e-06 1.81317e-06 ...\n",
      "    RERCLD               (lev, lat, lon) float32 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...\n",
      "    T                    (lev, lat, lon) float32 234.441 234.392 234.331 ...\n",
      "    TOT_CLD_VISTAU       (lev, lat, lon) float64 nan nan nan nan nan nan nan ...\n",
      "    TOT_ICLD_VISTAU      (lev, lat, lon) float64 nan nan nan nan nan nan nan ...\n",
      "    Z3                   (lev, lat, lon) float32 37596.6 37593.7 37591.5 ...\n",
      "    ACTNL                (lat, lon) float32 9.19856e+06 2.9595e+07 ...\n",
      "    ACTREL               (lat, lon) float32 1.13186 2.61381 6.92587 9.69679 ...\n",
      "    AEROD_v              (lat, lon) float64 nan nan nan nan nan nan nan nan ...\n",
      "    ASDIR                (lat, lon) float32 1.0 1.0 1.0 1.0 1.0 1.0 1.00186 ...\n",
      "    CDNUMC               (lat, lon) float32 2.50045e+09 2.9527e+10 ...\n",
      "    CLDBOT               (lat, lon) float32 25.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 ...\n",
      "    CLDHGH               (lat, lon) float32 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...\n",
      "    CLDLOW               (lat, lon) float32 0.361487 0.726613 0.714479 ...\n",
      "    CLDMED               (lat, lon) float32 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...\n",
      "    CLDTOP               (lat, lon) float32 22.0 30.0 30.0 30.0 30.0 30.0 ...\n",
      "    CLDTOT               (lat, lon) float32 0.361487 0.726613 0.714479 ...\n",
      "    FCTL                 (lat, lon) float32 0.161439 0.313795 0.714479 ...\n",
      "    FLNT                 (lat, lon) float32 271.24 274.347 269.555 268.142 ...\n",
      "    FLUT                 (lat, lon) float32 272.538 275.643 270.85 269.436 ...\n",
      "    LWCF                 (lat, lon) float32 14.7673 12.4005 18.1683 20.3177 ...\n",
      "    ODV_DST01            (lat, lon) float64 nan nan nan nan nan nan nan nan ...\n",
      "    ODV_DST02            (lat, lon) float64 nan nan nan nan nan nan nan nan ...\n",
      "    ODV_DST03            (lat, lon) float64 nan nan nan nan nan nan nan nan ...\n",
      "    ODV_DST04            (lat, lon) float64 nan nan nan nan nan nan nan nan ...\n",
      "    ODV_SSLT01           (lat, lon) float64 nan nan nan nan nan nan nan nan ...\n",
      "    ODV_SSLT02           (lat, lon) float64 nan nan nan nan nan nan nan nan ...\n",
      "    ODV_SSLT03           (lat, lon) float64 nan nan nan nan nan nan nan nan ...\n",
      "    ODV_SSLT04           (lat, lon) float64 nan nan nan nan nan nan nan nan ...\n",
      "    PS                   (lat, lon) float32 102307.0 102210.0 102121.0 ...\n",
      "    SWCF                 (lat, lon) float32 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ...\n",
      "    TAU_tot              (lat, lon) float64 nan nan nan nan nan nan nan nan ...\n",
      "    TGCLDCWP             (lat, lon) float32 0.0477388 0.0660905 0.0609831 ...\n",
      "    TGCLDLWP             (lat, lon) float32 0.0477388 0.0660905 0.0609831 ...\n",
      "Attributes:\n",
      "    Conventions: CF-1.0\n",
      "    source: CAM\n",
      "    case: aci_arg_comp_PD\n",
      "    title: UNSET\n",
      "    logname: darothen\n",
      "    host: ys0536\n",
      "    Version: $Name$\n",
      "    revision_Id: $Id$\n",
      "    initial_file: arg_comp_PD.cam.i.0011-01-01-00000.nc\n",
      "    topography_file: /glade/p/cesmdata/cseg/inputdata/atm/cam/topo/USGS-gtopo30_1.9x2.5_remap_c050602.nc\n",
      "    history: Wed Mar 30 19:08:08 2016: ncks --append aci_arg_comp_PD.cam.h1.0002-06-15-00000.nc medium_sample.nc\n",
      "    NCO: 4.4.7\n"
     ]
    }
   ],
   "source": [
    "test_ds = xr.open_dataset(\"medium_sample.nc\",)\n",
    "test_ds = test_ds.isel(time=0, lat=slice(-60, 60))\n",
    "\n",
    "print(test_ds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pressure-level Interpolation\n",
    "\n",
    "Select the 3D temperature field, and interpolate to specific pressure levels given the background hybrid pressure levels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "OrderedDict([('long_name', 'hybrid level at midpoints (1000*(A+B))'),\n",
       "             ('units', 'level'),\n",
       "             ('positive', 'down'),\n",
       "             ('standard_name', 'atmosphere_hybrid_sigma_pressure_coordinate'),\n",
       "             ('formula_terms', 'a: hyam b: hybm p0: P0 ps: PS')])"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "field = 'T'\n",
    "test_data = test_ds[field]\n",
    "\n",
    "# Compute the 3D pressure field given the hybrid coordinate system\n",
    "pres_levs = ma.hybrid_to_pressure(test_ds)\n",
    "\n",
    "# Select pressure levels for interpolation\n",
    "new_pres_levs = 100.*np.array([300., 500., 750., 850., 960.])\n",
    "# new_pres_levs = 100.*np.array(np.logspace(2, 3, 50))[-5:]\n",
    "\n",
    "# Compute interpolated data\n",
    "new_data_np = ma.interp_to_pres_levels(test_data, pres_levs, new_pres_levs,\n",
    "                                      'numpy')\n",
    "# new_data_sp = ma.interp_to_pres_levels(test_data, pres_levs, new_pres_levs,\n",
    "#                                       'scipy')\n",
    "\n",
    "# As a safety check, copy attributes into interpolated data\n",
    "nd = ma.copy_attrs(test_data, new_data_np)\n",
    "new_data_np.lev.attrs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Select a column to compare, and plot the original and interpolated data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATcAAAFTCAYAAABcR0t+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtcTPn/B/DXVNNFLZKwaEmWXIrKlohUyK2yKZFat3VJ\nF/fbsuKXjWWjRa0vya1cyrUlu265bVjVymKtpVB2o7tbl6n5/P5om21UTJqZM03v5+Mxj93OOTPz\nasx5d875nM/nw2OMMRBCiJJR4TqAIsrPz8fmzZuRn5/PdRQxipiLMkmGMklOWrmouNWgoKAAW7Zs\nQUFBAddRxChiLsokGcokOWnlouJGCFFKVNwIIUqJihshRCmpcR1AEZWWlgIAnj59Cj6fz3Ga/2Rl\nZYn+qyi5KJNkKJPknj59CuC//fBD8ehWkP/8/PPPOH36NJ49e4YbN25wHYeQRu2zzz5D69atMXTo\nUDg6Otb5+VTcavD48WMMHToU0dHRaNOmDddx6iQrKwvDhw/H9evXoa6uznUcQuosKysLEyZMwOnT\np9GhQ4cPfh06La2BqqoqAKBNmzZo3749x2nq5uDBg3B2dkanTp24jkJIvVTuhx+KGhSUzP79+zF+\n/HiuYxDCOSpuSuSvv/5CZmYm7OzsuI5CCOeouCmR/fv3Y+zYsfU+nCdEGVBxUxKMMezfvx/jxo3j\nOgohCoGKm5IYP3487t27h0OHDnEdhRCFQMVNCWRnZ+PgwYMAgB9++IHjNIQoBipuDZxAIMDYsWNh\naWkJTU1NzJo1i+tIhCgEus+tgVu0aBE0NTVx9uxZakggpAoqbg3Y3r178eOPP+LGjRtU2Ah5CxW3\nBiolJQXz5s1DQkICdHV1uY5DiMKha24NUHZ2NlxdXfHDDz+gZ8+eXMchRCFRcWtgysrK4OHhgfHj\nx8PNzY3rOIQoLCpuDczChQuhoaGB1atXcx2FEIVG19yqqBzP7c2bN1xHqRE1IJDG5JtvvkGTJk1o\nPDdpyszMhIODA86dO6cwQx6lpKTA0dER58+fh4mJCddxCJEZae1/dFraAFQ2IISHh1NhI0RCVNwU\nXNUGBHd3d67jENJgUHFTcIsWLYK6ujo1IBBSR9SgoMCioqIQFxeHX3/9lRoQCKkjKm4KKiUlBXPn\nzsX58+fRokULruMQ0uDQaakCysnJgaurK8LCwqgBgZAPRMVNwVQ2IIwbNw5jx47lOg4hDRYVNwWz\naNEiqKmp4ZtvvuE6CiENGl1zUyDR0dE4fvw49UAgRAqouCmI3377DXPmzKEGBEKkhE5LFUBOTg4+\n//xzakAgRIqouHGssgHBw8ODGhAIkSIqbhxbvHgx1NTUEBwczHUUQpQKXXOrQt5DHu3btw/Hjh2j\nBgRCalDfIY+ouFXh6OgIR0dHZGZm4vz58zJ9L29vb0RFRcHb25saEAipwbJly2jIo4Zo3759AIDY\n2FiOkxCinOjIjQNlZWVQU1ODiooKTaJMiIxQceNAcnIyunbtilu3bnEdhRClRaelHLhw4QLs7Oy4\njkGIUqPixoGEhAQqboTIGBU3ORMIBEhMTMTAgQO5jkKIUqPiJmdJSUkwMjKi2z8IkTEqbnJGp6SE\nyAcVNzmj4kaIfFBxk6PS0lJcu3YNAwYM4DoKIUqP7nOrQtZ9S3/99Vd07doVzZs3l8nrE6JMqG+p\nFMm6bymdkhIiOepb2oAkJCRg0KBBXMcgpFGg4iYnxcXFuHHjBl1vI0ROqLjJyfXr19G9e3c0bdqU\n6yiENApU3OSETkkJkS8qbnJCneUJkS8qbnJQVFSE5ORk2NjYcB2FkEaDipscTJkyBa9fv0ZgYCDX\nUQhpNKi4ycHhw4fBGEN4eDjXUQhpNKi4ycGwYcNoSHFC5Ix6KFQhq+5XHh4e0NHRQUhIiFRflxBl\nRt2vpEhW3a+Ki4uhpaUltdcjpDGg7lcNQHFxMTQ1NbmOQUijQsVNDoqKiujIjRA5U8jiduvWLbE+\nmC9evICfnx/69OkDe3t7HDp0SGz7kJAQWFtbw8rKCsHBwWCMyTvyO9GRGyHyp3DF7dChQ5g6dSrK\nyspEy5YvXw5tbW1cvXoVoaGhWL9+vWjOz6ioKFy6dAknTpxAfHw8kpOTERkZyVX8GtGRGyHyp1DF\nbevWrYiKioKPj49o2Zs3b3Du3DkEBASAz+fD1NQUTk5OOHbsGAAgLi4OEydOhJ6eHvT09DBjxgwc\nOXKEq1+hRnTkRoj8KVRxc3Nzw7Fjx9CzZ0/RskePHoHP56Ndu3aiZYaGhkhLSwMApKWloXPnzmLr\nHj16JLfMkigqKqLiRoicKVRxa9myZbVlRUVF0NDQEFumqamJ4uJi0fqqhUNTUxNCoRClpaWyDVsH\ndCsIIfKn8Pe5aWlpVStUxcXFaNKkCQDxQle5TlVVFerq6hK9fn5+PgoKCsSWZWVl1TO1ODotJaTu\nMjIyIBAIxJY1b94curq6Ej1f4Ytbhw4dIBAIkJWVhTZt2gAA0tPTYWRkBAAwMjJCeno6TE1NAVSc\nplauk0RUVBS2bNki/eBVUIMCIXU3adKkasv8/Pzg7+8v0fMVvrhpa2vD3t4eISEhCAoKwv3793Hi\nxAls374dAODs7IwdO3agb9++UFVVxbZt2zB69GiJX9/LywujRo0SW5aVlVXjB/uh6MiNkLrbtWuX\n6ICmUl1mjlP44gYAQUFBCAwMhK2tLbS1tbF48WKYmJgAADw9PZGbmws3NzcIBAK4uLjUqTDp6upW\nO8zl8/nSjI8//vgDLi4u8Pf3p/6lhEjIwMCgXt2veEzR7nhVAJmZmXBwcMC5c+fq9eFW4vP5KCsr\ng6amJoqKiqSQkBDlJa39T6FaS5WVh4cHeDweDXlEiBxRcZODPXv2oHXr1lTcCJEjKm5yoKKiguHD\nh+PUqVNcRyGk0aDiJicjRoxAfHw81zEIaTSouMnJkCFDcOXKFamP8ksIqRkVNzlp1qwZzM3NkZCQ\nwHUUQhoFKm5yNHLkSDo1JUROqLjJUeV1N7q1kBDZo+ImR927d4dQKMQff/zBdRRClB4VNzni8XjU\nakqInFBxkzMqboTIBxU3ObO3t0dSUhJevHjBdRRClBoVNznT1tZGixYt0LJlS8yfP5/rOIQoLSpu\nHHj69CkEAgHCw8O5jkKI0qLixgF/f3+oqqqiZcuW1GOBEBmh4saBDRs2oKSkBHZ2dnB2dqYx3giR\nASpuHFFVVcXOnTvRqlUrfP7552KT3BBC6o+KG4dUVVWxZ88eNG3aFG5ubigpKeE6EiFKg4obx9TU\n1BAdHQ11dXV4eHgo1HyrhDRkVNwUAJ/Px4EDByAUCjF+/PhqczUSQuqOipuCUFdXR2xsLIqLi+Hl\n5YWysjKuIxHSoFFxUyAaGho4fPgwCgoKMHHiRJSXl3MdiZAGi4qbgtHU1MSxY8eQlZWFqVOnQigU\nch2JkAaJipsC0tLSQlxcHB49eoTp06dTgSPkA1BxU1Da2to4ceIE/vzzT/j6+tIAl4TUUZ2KW3p6\nOpKSkpCSkoKMjAxZZSL/0tHRQXx8PG7evImAgAAqcITUgdr7NkhOTkZUVBR++eUXsWF6eDwemjVr\nhgEDBmD8+PEwNzeXadDG6qOPPsJPP/2EIUOGYN68ediwYQN4PB7XsQhReLUWt8ePH2PFihXIysqC\ng4MDvv/+exgZGaF58+YQCoXIz8/HvXv3kJSUhPnz56N9+/b4v//7PxgaGsozf6PQrFkz/Pzzzxg8\neDAsLS3x+++/w9fXFyEhIVxHI0RxsVp4eHiwS5cu1bZajFAoZGfPnmVjx46VaHtFl5GRwbp06cIy\nMjK4jiImJyeH8Xg8BoBpampyHYcQmZDW/lfrNbcDBw5gwIABEhVIHo8HBwcHHDx4UFo1l9RAT08P\nPj4+4PF4+PTTT+k+OELeoU4NCuXl5SgtLRU9Xr16hWvXrskqG6lBWFgYCgsLoaenhy+++IK6ahFS\ni/c2KADAr7/+ihUrVuDx48fVX0BNDb///rvUg5HaffTRR4iPj4ebmxvc3Nxw8OBBaGpqch2LEIUi\n0ZHbmjVr0KlTJ+zYsQNaWlrYvHkzli9fjqZNm2LdunWyzkhqoKWlhaNHj0JDQwOjRo3Cq1evuI5E\niEKRqLg9ePAA8+bNQ79+/dCjRw9oaGhgwoQJWLFiBSIjI2WdkdRCXV0d+/fvxyeffAJHR0cUFBRw\nHYkQhSFRcdPQ0IC6ujoAwNDQEPfu3QMAmJmZ4eHDh7JLR95LVVUVERER6NOnD+zs7JCdnc11JEIU\ngkTFzcLCAlu3bsWrV6/Qs2dPnDt3DuXl5UhJSUGTJk1knZG8h4qKCkJDQzFy5EgMHDgQT58+5ToS\nIZyTqLgtXrwYycnJiImJgbOzM169egULCwssXLgQXl5ess5IJMDj8bB69WpMnjwZAwYMQFpaGteR\nCOGURK2lnTp1ws8//4yioiJoamoiJiYGV65cQevWrdGrVy9ZZyR1sGjRInz00UcYOHAgTp8+je7d\nu3MdiRBOvLO4ZWRkICEhAXw+HwMGDED79u0BAE2aNMHQoUPlEpDUnY+PD3R0dODg4ICTJ09Sv1/S\nKNVa3C5evAhfX19oaGgAANauXYs1a9ZgxIgRcgtHPpy3tze0tbUxbNgwHD16FP379+c6EiFyVes1\nt02bNmHcuHG4fv06kpOT4eXlhfXr18szG6knV1dX7N27F6NHj8aZM2e4jkOIXNVa3B4+fIiJEydC\nTa3i4M7Hxwf//PMP8vLy5BaO1J+joyOOHDmCCRMm4Pjx41zHIURuai1uxcXFYrd56OjoQEtLC69f\nv5ZLMCI9AwYMQHx8PGbMmIERI0ZAU1MT8+fP5zoWITJFw4w3En369MHZs2dx6tQplJSUIDw8nOtI\nhMhUrcWNx+NBIBCIjQJS0zKaIb3h6NmzJyZPngwej4d+/fpxHYcQmaq1tZQxBjs7u2rLRo4cWW3b\nP/74Q/rJiExERkbi66+/hr29PTZu3Ii5c+dyHYkQmai1uO3Zs0eeOYgcGRoa4uLFi7C3t0dxcTGW\nLl3KdSRCpK7W4mZpaSnPHETOPvnkE1y6dAkODg4oKSlBYGAgTTxDlIpE3a9KS0sRGxuLu3fvori4\nuNp6mqikYWrbti0uXLiAwYMHo6SkBMHBwVTgiNKQqLgFBgbi1KlT6N+/P5o2bSrrTESOWrdujYSE\nBAwZMgQlJSUICQmhAkeUgkTF7aeffkJYWBh14VFSLVu2xPnz5+Ho6Ag/Pz9s3rwZKip0lxBp2CT6\nBmtra+Pjjz+WdRbCIV1dXZw5cwY3b97EjBkzaGYt0uBJVNy++OILbNiwgcbpV3KVkz//9ddfmDx5\nMsrKyriORMgHq/W01MbGRuznnJwcWFlZoVmzZtVOWa5cuSKbdETudHR0EB8fj9GjR8PLywt79+4F\nn8/nOhYhdVZrcZs3b57cLywnJSVh3bp1SEtLQ4sWLTB16lR4eHjgxYsX+Oqrr3Dt2jU0bdoUs2bN\ngpubm+h5ISEhOHToEIRCIVxcXLB06VK6KF4PTZo0QVxcHMaMGQMPDw8cOHBANIcGIQ1GbVPRC4XC\nOk9f/yHPqVRYWMgsLS3ZyZMnGWOM3blzh1laWrLExETm7+/PFi1axEpLS1lqaiqztLRkqampjDHG\n9u7dy5ydnVlOTg7Lyclhrq6uLCIi4oNzMMZYRkYG69KlC8vIyKjX6zR0JSUlbPTo0WzkyJGsqKiI\n6zikkZDW/lfrNTc3NzecPXtWogIpFAoRHx+PMWPGfHCR/fvvvzFo0CDRYJjdu3eHlZUVUlJScP78\neQQEBIDP58PU1BROTk44duwYACAuLg4TJ06Enp4e9PT0MGPGDBw5cuSDc5D/qKurIyYmBtra2nB2\ndsabN2+4jkSIxGotbqGhoYiJiYGdnR1Wr16NhIQEZGRk4NWrV3j58iWePHmC06dP45tvvsHAgQNx\n5MgRhIaGfnAQY2NjfPvtt6KfCwsLkZSUBKBiVvt27dqJ1hkaGoomQElLS0Pnzp3F1j169OiDcxBx\nfD4f0dHRaNOmDbp27UrDJZEGo9ZrbgYGBti2bRtu3bqF6OhofPXVV8jPzxe7lqWnp4eBAwciPDwc\npqamUgv18uVL+Pj4wMTEBFZWVtX6uWpqaop6SlROWlN1nVAoRGlpqUTXifLz86tNZpyVlSWF30J5\nqKmpYefOndDQ0EB5eTnCw8OpVwqRuYyMDAgEArFlzZs3h66urkTPf+9NvKampqLC9fTpU+Tm5kJF\nRQX6+vpo3br1B0R+t4yMDPj4+KBDhw7YuHEjHjx4UG1YpaoDaVYtdJXrVFVVJb4AHhUVhS1btkjv\nF1BSqqqqCAgIwObNm8Hn83H79m307NmT61hEiU2aNKnaMj8/P/j7+0v0fIl6KFRq166d2OmhtN25\ncwfTpk2Di4sLFi9eDADo0KEDBAIBsrKy0KZNGwBAeno6jIyMAABGRkZIT08XFeC0tDTROkl4eXlh\n1KhRYsuysrJq/GAbuw0bNmDDhg2IioqCvb09oqOjMWTIEK5jESW1a9cu0T5fqXnz5pK/gJQaOOot\nOzubWVtbs+3bt1db5+/vzxYsWMCKiopYamoqs7KyYrdu3WKMVbSWOjk5saysLJadnc1cXV1ZZGRk\nvbJQa+n7Xbp0ibVu3Zpt27aN6yhEyUhr/6vTkZssHT58GPn5+QgPD0dYWBiAitGAv/jiC6xevRor\nVqyAra0ttLW1sXjxYpiYmAAAPD09kZubCzc3NwgEAri4uNBRlxwMGDAAly9fxogRI/DXX39h7dq1\n1B+VKBQeY4xxHULRZGZmwsHBAefOnRNNRE1qlpubi88//xytWrXCnj17xCYVIuRDSGv/q9Of2mfP\nnuHatWsoLi5GTk7OB78pUR56eno4c+YMNDU1YWdnRy3NRGFIVNzevHmDuXPnwtbWFlOmTEF2djZW\nrFgBT09PmseUQENDA3v37sWIESNgbW2NO3fucB2JEMmK2/r165GVlYVTp05BQ0MDADB//nzR6K2E\n8Hg8BAYGIigoCHZ2djTDPeGcRMXt3LlzWLp0KQwNDUXLjIyMsGrVKly+fFlm4UjD4+XlhUOHDsHb\n2xvbt2/nOg5pxCRqLX316hV0dHSqLVdRUaExv0g1AwcOxOXLlzFy5Eg8ePAAa9asoZZUIncSfeNs\nbGywdetWsdFZ8/PzsX79ehp6nNTo008/xdWrV5GYmIixY8eiqKiI60ikkZGouC1fvhyPHj2CtbU1\niouL8eWXX8LOzg6FhYVYtmyZrDOSBkpPTw9nz56FpqYmBg0ahGfPnnEdiTQiEp2WNm3aFDExMbh6\n9SrS0tJQVlYGIyMj9O/fnwaFJO9U2ZK6atUq9O3bFydPnkT37t25jkUaAYmKm5OTEzZt2gRra2tY\nW1vLOhNRMjweDytXroSRkREGDRqEffv2YfDgwVzHIkpOotNSgUBAR2ik3ry9vXHo0CG4uLiAz+dj\n3rx5XEciSkyiI7dRo0Zh8uTJGDlyJAwMDMTGTwMADw8PmYQjymfgwIEoKytDWVkZvv/+eyxfvhwt\nWrTgOhZRQhIVt/j4eGhpaeH8+fPV1vF4PCpupE78/PwQHh6OHj16wMzMDAcOHKDLHUTqJCpuNRU1\nQj5USEiIaCTfuLg4jB49GgsWLMD8+fPpfjgiNRIVtxs3brxz/WeffSaVMKTxcXZ2Rq9evTBu3Dhc\nvHgRu3btQsuWLbmORZSARMXN29u7xuV8Ph/NmjWjSZlJvXTo0AGXLl3CsmXLYG5ujv3799PN4aTe\nJCput27dEvu5rKwMT548QXBwMF1vI1LB5/Oxbt062NraYsyYMZgzZw4WLVpEp6nkg0n0zVFXVxd7\nNGnSBMbGxvjqq69oFiQiVSNHjsSNGzfw448/YuTIkcjOzuY6Emmg6vVnsaioCPn5+dLKQgiAimkl\nL1y4gF69esHc3JxGniEfRKLT0g0bNlRb9urVK5w+fRoDBgyQeihC+Hw+1q5di4EDB8Ld3R0BAQFY\nsmQJnaYSiUlU3H777Texn3k8Hvh8Ptzc3DBlyhSZBCMEAEaMGIGkpCSMHz8eFy9exN69e9GqVSuu\nY5EGQKLitnfvXlnnIKRW7du3R0JCAlasWAFzc3NER0fD1taW61hEwUl8jH/+/HnRxd2DBw9i4sSJ\nWL9+fbXZ4AmRBTU1NQQHB2PHjh0YN24cVq9eLTa+ICFvk6i4hYWFYf78+Xj69Cl+/fVXrFq1Ch07\ndsSFCxewZs0aWWckRMTR0RFJSUk4c+YMhg0bRmPEkVpJVNwOHTqE0NBQ9O7dG3FxcbCwsMCqVauw\nZs0anDp1StYZCRHTrl07nDt3Dn379oWFhQUSEhK4jkQUkETFLS8vD126dAEAXLhwAYMGDQIANGvW\njE5LCSfU1NQQFBSEnTt3wsnJCXw+H/Pnz+c6FlEgEhU3IyMjxMbGIjo6Gjk5OXBwcEBxcTG2bduG\nHj16yDojIbUaMmSIaAilTZs2gTHGdSSiICQqbkuXLkV0dDSCgoLg5eWFjh074ttvv8WlS5fw1Vdf\nyTojIe/k6+sLTU1NtGzZElOmTIFAIOA6ElEAEhW3zz77DFevXsX169exfPlyABVfqISEBHTr1k1q\nYeLj4zFixAiYmZnByckJZ8+eBQC8ePECfn5+6NOnD+zt7XHo0CGx54WEhMDa2hpWVlYIDg6mv96N\nTEhICIqKivDgwQM8e/YMn3/+Od68eSO+UWxsxYM0HkxC9+7dYwUFBYwxxi5cuMCWLl3KoqOjJX36\ne6Wnp7PevXuzmzdvMsYYS0xMZD179mT5+fnM39+fLVq0iJWWlrLU1FRmaWnJUlNTGWOM7d27lzk7\nO7OcnByWk5PDXF1dWURERL2yZGRksC5durCMjIx6/15EvkpLS5mXlxeztrZmubm5FQufP2dMX7/i\n8fw5twHJe0lr/5PoyO3gwYP4/PPP8eeff+Lu3bvw9fXF8+fPER4ejtDQUKkU2Y4dOyIxMRG9evVC\nWVkZsrOzoaOjAzU1NZw7dw4BAQHg8/kwNTWFk5MTjh07BqBisMOJEydCT08Penp6mDFjBo4cOSKV\nTKTh4fP52L17N/r37w8bGxtkZGQAs2YB2dkVD19friMSOZGouEVGRmLNmjWwtLTEkSNHYGxsjIiI\nCGzYsAGHDx+WWhgtLS1kZmaiV69eWLJkCebOnYuMjAzw+Xy0a9dOtJ2hoSHS0tIAAGlpaejcubPY\nukePHkktE2l4VFRUsH79ekyZMgXB5uZA1csYdHraaEhU3P755x9YWloCqLgVxM7ODgDQtm1bvHr1\nSqqB2rZti1u3bmHnzp1Ys2YNzp8/Dw0NDbFtNDU1UVxcDKBiZJKqE9ZoampCKBTSLSoECyZOxMaS\nkuorfH0rjuKIUpOob6mBgQEuXbqEVq1aITMzE/b29gCAw4cPo1OnTlINVDnqg5WVFRwdHXH79u1q\nhaq4uBhNmjQBIF7oKtepqqpCXV1dovfLz89HQUGB2LKsrKz6/ApEUcyaBc2XL6svrzw9jYmRfyYi\nsYyMjGot382bN4eurq5Ez5foyG327NkICgrCrFmz4ODggG7duiE4OBi7du3CggUL6p66BhcvXsTk\nyZPFlgkEAnTo0AECgUCs4KSnp8PIyAhAxT146enponVpaWmidZKIiorCsGHDxB6TJk0CUHGKy+Px\nxB4rV66s8XVWrlxZbVvantvtY99qVa/qzp07Cp+/sW8/adKkavtmVFRUja9XI0lbHnJzc9ndu3dF\nPz98+JDl5OTUqzWjquzsbPbZZ5+x48ePM6FQyC5cuMD69OnD0tLSmL+/P1uwYAErKipiqampzMrK\nit26dYsxVtFa6uTkxLKyslh2djZzdXVlkZGREr9vXl4eS0tLE3skJiZSa6kyqGwlBcQer7S1mfDZ\nM67TkVpUtpYmJiZW2zfz8vIkfh2Ji1tJSQmLi4tj33//PcvPz2fXrl1j2dnZHxS+NklJSczV1ZVZ\nWFiwMWPGsF9//ZUxxlhBQQGbPXs2s7S0ZHZ2duzIkSOi55SXl7PQ0FBmY2PDrKysWHBwMBMKhfXK\nQbeCKJGYmGrFbW779mzu3LmsvLyc63SkBtLa/yQqbk+ePGF2dnZs4MCBrHv37uzJkyfM39+fWVpa\nsjt37tQrgCKi4qZk3Nz+K27u7iwvL4/179+feXp6spKSEq7TkbfI9T63b775BjY2NkhISBBdqN+w\nYQMGDRpEQx4RxRceDujrVzzCwqCrq4szZ87g1atXcHZ2lnqLP1EMEhW35ORkTJo0SWz8ejU1Nfj4\n+OD27dsyC0eIVPxb1BAWVvH/qLin8vDhw2jXrh3s7e1pli0lJPHUfi9evKi2PDMzE9ra2lIPRYjU\nubtXPKpQU1NDREQEeDweWrdujWnTpnEUjsiCRMXN2dkZQUFBoqO0wsJCXLx4EStWrMCoUaNkGpAQ\nWeLxeEhNTQVjDDt27MD9+/e5jkSkRKLiNn/+fFhZWcHT0xNFRUVwc3ODr68vHBwcMG/ePFlnJESm\nKodMGjJkCAYNGoTU1FSuIxEpkKiHwu+//445c+Zg9uzZePLkCcrLy/HJJ5+IegkQ0pCFhIQgJCQE\nABAbG4uhQ4fi6NGj6NevH8fJSH1IdOTm4+ODhw8fQkNDA59++imMjY2psBGl5O7ujt27d2P06NE4\nc+YM13FIPUhU3Nq3by/WxYkQZTZs2DAcOXIEEyZMwNGjR7mOQz6QRKelRkZGWLBgAbZu3QoDAwOx\nUTgAiA7pCVEWNjY2+OmnnzBy5Ei8fPkSX3zxBdeRSB1JVNxUVFTg4uIi6yyEKBRzc3MkJCRg6NCh\nKCwshL+/P9eRSB1IVNyoFwJprIyNjXHp0iUMGTIEhYWFWLZsGXg8HtexiARqLW7l5eXYtm0bTp8+\nDT6fj8GDB2Py5Mng8/nyzEcI5zp27IjLly9j6NChKCgowPr166nANQC1Nihs3LgR27dvh6mpKUxM\nTLBt2zb83//9nzyzEaIw2rRpgwsXLuDKlSuYPn06ysvLuY5E3qPW4nby5EmsX78eq1atwtdff43N\nmzcjLi4na1TzAAAgAElEQVSO/lFJo9WiRQucPXsWaWlp8PT0pKHsFVytxS07Oxs9e/YU/WxlZYWy\nsjLk5OTIJRghikhHRwcnT55ESUkJRo8eXX1+VKIwai1uZWVlUFP775KciooK1NXV6a8VafQ0NTUR\nGxsLPT09DB8+vMZBJQj3JLqJlxAirnJ+1J49e8Le3p7OaBTQO28FOX78uNiQRkKhECdOnECLFi3E\ntvPw8JBqqJycHDg7O2PNmjWwtbXFixcv8NVXX+HatWto2rQpZs2aBTc3N9H2ISEhOHToEIRCIVxc\nXLB06VJqzSIyp6Kigi1btqBv375o1aoVvvzyS2zbto3rWORftRa3tm3bVptpRk9PD7FvTWjL4/Gk\nXtyWLVuGwsJC0c/Lly+HtrY2rl69ij/++APTpk1Dly5dYGpqiqioKFy6dAknTpwAAEyfPh2RkZGY\nOnWqVDMRUpOqQyZFRERg4cKF+PTTT7mORfCO4nb+/Hl55hA5cOAAtLW10aZNGwDAmzdvcO7cOdH9\ndqampnBycsKxY8dgamqKuLg4TJw4EXp6egCAGTNm4Pvvv6fiRuTG19cX4eHhGDBgAGxtbREfH4/e\nvXtzHavRU6hrbunp6di5cydWrlwJxhgA4PHjx+Dz+WjXrp1oO0NDQ6SlpQGomKe0c+fOYusePXok\n19ykcQsJCUFRURFOnz6NTZs2wdHREVeuXOE6VqMnUfcreSgvL8fixYvx9ddfo2nTpqLlb968gYaG\nhti2VWeZLyoqEuvIr6mpCaFQiNLSUolmnacZ54k0ubm5oVmzZnB1dcWuXbswYsQIriM1WPWdcV5h\niltYWBi6desGGxsbseVaWlrVbj8pLi4WjSdXtdBVrlNVVZWosAEVM85v2bKlnukJ+c+QIUMQFxcH\nFxcXhIaGYvz48VxHapAmTZpUbZmfn5/EAxgoTHE7deoUcnJycOrUKQDAy5cvMXfuXHz55ZcQCATI\nysoSXYdLT0+HkZERgIrhmNLT02Fqagqg4jS1cp0kvLy8qs0DkZWVVeMHS4ik+vbti3PnzmHYsGEo\nKCiAj48P15EanF27don2+UrNmzeX+PkKVdyqsre3R2BgIGxtbXHv3j2EhIQgKCgI9+/fx4kTJ7B9\n+3YAFZPX7NixA3379oWqqiq2bduG0aNHS/y+urq61Q5zaXAAIg09e/YUjSiSl5eHr776im5RqgMD\nAwO0b9/+g5+vMMXtbVW/BEFBQaJCp62tjcWLF8PExAQA4OnpidzcXLi5uUEgEMDFxYWOuojC6NSp\nE65cuYKhQ4ciLy8P3333HRU4OeGxymZJIpKZmQkHBwecO3euXn85CKmUn5+PkSNHwtjYGNu2bRPr\n2kjESWv/U6hbQQhRVrq6ujhz5gyePn0Kd3d3sUYwIhtU3AiRE21tbcTFxYHP54vmZiCyQ8WNEDnS\n0NDA/v37YWRkBAcHB+Tm5nIdSWlRcSNEzlRVVfG///0P9vb2GDhwIJ4+fcp1JKVEVzUJ4QCPx8Pa\ntWuhq6sLGxsbnD59mjrcSxkduRHCocWLF+Orr76CmZkZNDQ0MH/+fK4jKQ0qboRwbNq0aSgpKUFp\naSnCw8O5jqM0qLgRogACAgKgpqaGZs2a4fXr11zHUQpU3AhRACEhISgtLcWwYcMwfvx4lJWVcR2p\nwaPiRoiC4PF42LZtG4qKihAQEADqPFQ/ClXcIiMj0bNnT5ibm8PMzAzm5uZITk7Gixcv4Ovriz59\n+sDe3h6HDh0Se15ISAisra1hZWWF4OBg+lKQBktdXR2HDx/GL7/8gnXr1v23Ija24kEkplC3gty9\nexcLFiyo1vE9ICAAOjo6NIcCaRSaNm2K+Ph4WFtbw8DAAJ5DhgC+vhUrBw0C9PU5zddQKNSR2x9/\n/IGuXbuKLaucQyEgIKDaHAoAxOZQ0NPTw4wZM3DkyBEu4hMiNe3atcPJkycxZ84cPHd3B7KzKx6V\nRY68l8IUt+LiYqSnp2PPnj2wsbHByJEjcfjwYZpDgTRaJiYmOOfjg1YXL/63kE5PJaYwp6U5OTmw\nsLCAp6cnrK2tcfPmTfj4+GDy5MkynUOBEIWVnQ2TH36ovtzXl05PJaAwxa19+/bYu3ev6Oc+ffrA\nxcUFSUlJMp1DgSaIIQpr1qyKU9G3VZ6exsTIP5Mc1XeCGIU5Lb1792612bpLSkrQtm1b0RwKlWqa\nQ6FSXedQiIqKwrBhw8QelQ0ahoaG4PF4Yo+VK1fW+DorV66sti1tT9vXZ/vYt+4KqOrOnTsKn7++\n20+aNKnavvn2RPHvxBREeno6MzU1ZT///DMTCoUsMTGRmZubs7t37zJ/f3+2YMECVlRUxFJTU5mV\nlRW7desWY4yxvXv3MicnJ5aVlcWys7OZq6sri4yMlPh98/LyWFpamtgjMTGRdenShWVkZMjq1yXk\n/Z4/Z0xfnzFA7FGoqcmEz55xnU5mMjIyWJcuXVhiYmK1fTMvL0/i11GY4sYYYwkJCczJyYn17t2b\nDRs2jJ0+fZoxxlhBQQGbPXs2s7S0ZHZ2duzIkSOi55SXl7PQ0FBmY2PDrKysWHBwMBMKhfXKUfnh\nUnEjnIuJqVbcFnfqxAIDA7lOJjPS2v8UqrgpCipuRKG4uf1X3NzdWVZWFjM0NGQ7duzgOplMSGv/\nU5hrboSQWoSHV7SM6usDYWFo3bo1Tp06BX9/f6irq9MwSbWg4kaIovu3qCEsTHT7R9euXSEQCCAQ\nCGiYpFpQcSOkIXB3r3hU4e/vD1VVVRgYGFB/6hpQcSOkgQoJCcGrV6+gpaWFnTt3ch1H4VBxI6QB\n09TUxL59+7B48WL89ddfXMdRKFTcCGngevTogcDAQHh6ela7o78xo+JGiBLw9fVFq1atau0R0Bgp\nVHF79uwZZs6cCQsLCwwaNEjU1/TFixfw8/OjwSoJqQWPx0NkZCR27tyJS5cucR1HIShUcZs1axY6\nd+6MGzduYMeOHdiyZQtu3ryJ5cuXQ1tbG1evXkVoaCjWr1+PW7duAYDYYJXx8fFITk5GZGQkx78J\nIfLXunVrREREwNvbu9pgEI2RwhS31NRUZGdnY/78+VBRUYGRkREOHjyIVq1a0WCVhEhoxIgRcHZ2\nxsyZMxv9GYzCFLc7d+6gc+fOWLduHWxsbDBs2DDcvHkThYWFNFglIXWwbt063L59W2wIscZIYcZz\nKywsxPXr12FtbY0LFy7g999/x7Rp07B161aZDlZJ47kRZaOlpYV9+/bBwcEB/fv3r9MQYIqkvuO5\nKUxxU1dXR/PmzTFt2jQAgJmZGYYMGYLNmzfLdLDKqKgobNmyRUq/BSGKwdTUFMuWLYOXlxcuX74M\nNTWF2dUl9vZEUQDg5+cHf39/iZ6vML+xoaEhysrKwBgDj8cDAAiFQnTv3h3JycnIyspCmzZtANQ8\nWKWpqSmAug9W6eXlhVGjRokty8rKqvGDJaQhCQgIwKlTp7B69eoGeYvIrl27RPt8pebNm0v+AlIY\noUQqiouLma2tLdu0aRMrKytjycnJzNzcnKWmpsp0sMqa0JBHRFn8/fffrEmTJozP57N58+ZxHUci\nSjfkkYaGBvbu3YvU1FT069cPCxcuxNdffw1TU1MEBQVBIBDA1tYWc+bMweLFi2FiYgIA8PT0hIOD\nA9zc3DBq1Cj06dOHjroI+dfHH38sGj0kLCyM6zhyxWOskbcX1yAzMxMODg44d+4c2rdvz3UcQupl\n/vz52LRpE/T19fHo0SOFnxVOWvufwhy5EUJkIyQkBCUlJbC0tMTcuXO5jiM3VNwIaQRUVFSwe/du\nnD17Frt37+Y6jlwoTGspIUS2mjVrhqNHj8LW1hampqYwMzPjOpJM0ZEbIY1I9+7dER4eDldXV+Tm\n5nIdR6YUprj9+OOPMDMzg7m5OczNzWFmZoZu3bphxYoVePHiBXx9fWlUEEKkwN3dHWPGjMGECRNQ\nXl7OdRyZUZji5uTkhN9++w0pKSlISUlBeHg49PX14evri+XLl0NHR4dGBSFEStauXYuSkpIGeXOv\npBSmuFX1+vVrLFmyBCtXrsRHH31Eo4IQImVqamo4ePAgdu/ejbi4OK7jyIRCFreIiAh07doV9vb2\nePz4MY0KQogMtGrVCrGxsfjyyy9x//59ruNIncIVtzdv3iA6Ohp+fn6inz90VBBCyLtZWVkhKCgI\n/fr1g6amplJN8Kxwt4KcPXsW7dq1E3WE19LSkumoIDTkEWnspk+fDl9fX5SXlyM8PBwhISFcRwKg\nREMeVUpISMDw4cNFP3fo0AECgUBmo4LQkEeksePxePDz88PmzZvRq1ev/1bExlb8963JoOWlvkMe\nKcyoIJXs7OzY9evXxZbJclSQvLw8lpaWJvZITEykUUFIo/P8+XNmaGjIoqKiGHv+nDF9/YrH8+dy\nzVE5KkhiYmK1fTMvL0/i11GoIzehUIisrCzo6+uLLQ8KCkJgYCBsbW2hra1dbVSQ3NxcuLm5QSAQ\nwMXFpU6jgujq6lY7zOXz+fX+XQhpaPT19REXFwd7e3sM69kTetnZFSt8fYGYGLnnMTAwqFfHeRoV\npAY0KghpzJKXLoXF2rXiC2Ni5HZ6SqOCEEKkLzsbFjt2VF/u6wtUHsk1EFTcCCH/mTWr5iKWnV1R\n4BoQhSpuKSkpGDNmDCwsLDB8+HCcOHECAHczzhsaGoLH44k9auuusnLlymrb0va0fUPdviYxsbEK\nlee9ZNXiUVfl5eXM2tqanT59mjHG2I0bN1iPHj3Y06dPmb+/P1u0aBErLS1lqampzNLSkqWmpjLG\nKlpLnZ2dWU5ODsvJyWGurq4sIiKiXlloDgXSaFW2kgJij2cA+3rmTLlEULo5FF68eIH8/HzRTXs8\nHg98Ph8qKirUt5QQedHXB2qYa8EXwLcNbEAKhSluzZs3x/jx4zFv3jz06NED3t7eWLFiBfLz86lv\nKSHy5O4OuLmJfrzZpQvi1NWhqqqKpKQkDoPVjcIUN8YYNDU1sXnzZqSmpuKHH37AN998g1evXlHf\nUkLkLTy84ihOXx+9r1xBSUkJ9u3bBycnpwbTyV5hbuI9ffo0fv/9dyxatAgAYGtri0GDBsl8xnnq\nW0pIDaqenv57U/3o0aORm5sLR0dHXLlyRexsShaUpm/pP//8U62IqampoUePHkhJSaG+pYTIWw03\n7U6dOhU5OTlwdHTEpUuX0KJFi4oVMuiHqjR9S//8809mYmLCjhw5whhj7Pr168zCwoLdvn2b+pYS\nokCEQiGbN28e69evH3v9+rXU+6FKq2+pwhQ3xhhLSEhgLi4uzMLCgo0aNYqdPXuWMcZYQUEBmz17\nNrO0tGR2dnaiAshYxS0koaGhzMbGhllZWbHg4GAmFArrlYNuBSHk3crLy5m3tzczNDRkh1RU/rtt\nxN293q8trf2P+pbWgPqWEvJ+AoEAXhoaOPh2CalnP1TqW0oI4RS/oAARVe5UEFGQfqgKVdyuXr2K\nzz//HBYWFhg3bpxohiuuul8RQt5h1ix8VFRUfbmi9EOt9wmylGRmZrLevXuz2NhYVl5ezi5cuMCs\nrKxYTk4OZ92v1NTUGACxR2BgYI3PCQwMrLYtbU/bK/X2bm7VumlVPg7W4/Wldc1NYYrbvn37mIeH\nh9iy6dOns71797IePXqwzMxM0fKgoCC2atUqxhhj7u7u7PDhw6J1P//8MxsxYkS9slCDAiESeEc/\n1PYaGh/8skrXt1QoFIr1NAAq+pdeuXIFampq1P2KEEVTSz9UPwBdbWzkn+ctClPcbGxskJqaitOn\nT6OsrAyXLl3C1atXUVRURN2vCFFUb/VDhbs71j58iGvXrkFdXZ3TqQIVpodChw4dEBoaig0bNiAw\nMBD9+/fH8OHDkZ2dTd2vCFFk4eHAxYsV/x8Whk76+hAIBBAIBNiyZcsHTxWoNN2vXr9+jY8//hjH\njx8XLfPw8MCECRNw/fp16n5FiKKqoR+qn58ftmzZAlVVVZw5cwZDhgz5b3sJu2opTferytbSO3fu\nsJKSEhYVFcUGDhzIioqKqPsVIQ3U5cuXmaamJuPz+WzevHkSddVSyu5XcXFxzN7enllYWDBvb2/2\n4MEDxhh1vyKkIVNXV2cAGJ/PF799pJauWtT9Soao+xUh0jN//nyEhYXh87Iy7C8vF19ZQ1etBt/9\n6tatWxgwYIDo5/r0Qjhx4gQGDx4MMzMzzJw5E7m5uXL7PQgh7xYSEoLijAx8/3ZhA2TaVYuT4nbo\n0CFMnToVZWVlomXLly+HtrY2rl69itDQUKxfv17U/SoqKgqXLl3CiRMnEB8fj+TkZET+O577vXv3\nsHLlSmzcuBHXr19Hy5YtsXTpUi5+LUJIbWbNQqualsuwq5bci9vWrVsRFRUFHx8f0bI3b97UeRKY\no0ePAvjvqM3ExATq6upYsGABLl++jLy8vHpnpan9aHvaXvZT9cXExirH1H7Z2dmMsYrBKPv27csY\nY+zOnTusV69eYttFRUWxiRMnMsYYs7CwEPUlZYyx+/fvs+7duzPGGPPx8WHbtm0Te66VlRW7cePG\nB2ekBgVCpKyWrlrPeTzWtrIl9V8NtvtVy5Ytqy2rTy+EoqIiaGlpiT1XS0tL7MZeQgjHaumqFaCi\ngr8FAmzatEnqvYoU4iZeLS2tD+6F8PY6oKIYVj73faiHAiFyUtlVq7Kx0N0dbQ0MoBkeDg0NDTRp\n0gSzZ8/G3LlzAShJD4UOHTpAIBB8UC+EynWV8vLy8OLFC4l7KVAPBULk6K2uWiH6+ggJCYGmpibK\ny8sRHh4uKm717aGgEMVNW1sb9vb2CAkJQVBQEO7fv48TJ05g+/btAABnZ2fs2LEDffv2haqqKrZt\n24bRo0cDAEaNGgVvb2+MGTMGPXr0wIYNGzBw4EA0a9ZMovf28vLCqFGjxJZlZWXV+MESQuqphq5a\nAODr64vw8HDMmjVLtGzXrl2ig51KzZs3l/itFKK4AUBQUBACAwNha2sLbW1tLF68GCYmJgAAT09P\n5Obmws3NDQKBAC4uLqLiY2xsjKCgICxduhS5ubno06cPgoODJX5fXV3daoe5KioVlyLp9JQQGbC2\nrvhvZqZo0dy5c0VHbJX7Xdu2bdGhQ4cPfhvqoVCDhIQEzJw5k+sYhDRqW7duhZ2d3Qc/X2GO3BRJ\n27ZtAQA7duxAx44duQ1TRUZGBiZNmoRdu3bBwMCA6zgAKJOkKJPkHj16hKlTp4r2ww9Fxa0GlWPB\ntWvXTqH6lla2HLVp00ZhclEmyVAmyVXmknRMxtoozEi8hBAiTVTcCCFKiYobIUQpqa6Uak9V5aGp\nqQlLS8tqXbu4poi5KJNkKJPkpJGLbgUhhCglOi0lhCglKm6EEKVExY0QopSouBFClBIVN0KIUqLi\nRghRSlTcCCFKqVEXt7fnTn3+/DlmzpwJS0tLDBgwABs3bhTb/l1zp9ZXUlISxo4diz59+mDo0KE4\nePAggPrN5yqrTM+ePYOvry+srKxgY2OD1atXiw0HzUWmSowxeHt7Y926dWLLZZnpXbkEAgGCgoLQ\nt29f9O3bF8uXL+f8s+Lyex4fH48RI0bAzMwMTk5OOHv2LAAZfc/rNb1MAxYbG8v69OkjmoGLMcb8\n/f3ZmjVrmFAoZFlZWczBwYEdO3aMMcbY3r17mbOzM8vJyWE5OTnM1dWVRURESCVLYWEhs7S0ZCdP\nnmSMVcwGZmlpyRITE5m/vz9btGgRKy0tZampqczS0lI0ExhXmby8vFhQUBArLS1lOTk5bOzYsSw0\nNJTTTJW2b9/Ounfvzr799lvRMllmel+uNWvWsIkTJ7IXL16wwsJC5uHhwf73v//JPNf7vlNcfM/T\n09NZ79692c2bNxljjCUmJrKePXuy/Px8mXzPG+WRW01zpwIV8zaUlZWhrKwMjDGoqqqKun/UNHfq\nkSNHpJLn77//xqBBgzBixAgAQPfu3WFlZYWUlBScP3++TvO5yiOTtrY2fHx8wOfzoaenBycnJ/z2\n22+cZap873v37uHo0aMYPHiw2PNkmelduZKTkxETE4Ovv/4aH330EZo2bYrNmzfDyclJ5rne9e/3\n6NEjTr7nHTt2RGJiInr16oWysjJkZ2dDR0cHampqdZ63WJJMjbK4ubm54dixY+jZs6fY8i+//BIx\nMTEwMzODnZ0dzM3NMXToUAAVk9J07txZtK2hoSEePXoklTzGxsb49ttvRT8XFhYiKSkJAKCmpoZ2\n7dqJvW9aWhpnmbp3746tW7dCT09PtC4hIQHdunXjLFO3bt1QWlqKJUuWYPXq1dVmPpNlpnflatas\nGYRCIVJTU+Ho6AhbW1vs3LkTrVq1knmud31WU6dO5eR7DlTMdJeZmYlevXphyZIlmDt3LjIyMsDn\n86X+PW+Uxa2muVOBius1M2fOREpKCk6cOIGkpCTExMQAePfcqdL08uVL+Pj4wMTEBFZWVh88n6u0\nM82cORMmJibVhn1evXo10tPTMX36dE4zVU4MZGZmVm1beWWqzFX571dZdC9cuIDDhw8jJiYGv/zy\ni2jiIy4+K3t7e86/523btsWtW7ewc+dOrFmzBufPn5fJ97xRFreaPH/+HCtXrsS0adOgrq4OIyMj\nTJ8+XXQR9l1zp0pLRkYGxo8fD11dXWzevBlNmjT54PlcpZ2pRYsW2Lx5s2h5SUkJAgIC8MsvvyAq\nKko0yQ4Xma5evYpr164hICCgxu3lkalqrsp/Pw0NDTDGMGfOHOjo6KB169aYPHmy6CI6F5+VInzP\nVVRUoKqqCisrKzg6OuL27dsy+Z5TcftXTk6O6DpEJRUVFaipVYzE/vb8qFXnTpWGO3fuwMPDAwMG\nDEBYWBjU1dXF5nOtVNN8rvLMBFSc4nh5eeHly5eIiYkRG+uei0ynTp1CRkYG+vXrB0tLS/z444+I\njo4WTfIj60y15erYsSNUVFTEdtzK61zyyFVTJi6/5xcvXsTkyZPFlgkEAtl9z6XSDNJAXb9+XdRa\nWlJSwuzt7dmKFStYSUkJy8jIYCNHjmS7d+9mjFW02Dg5ObGsrCyWnZ3NXF1dWWRkpFRyZGdnM2tr\na7Z9+/Zq6/z9/dmCBQtYUVERS01NZVZWVuzWrVucZRIKhczLy4v5+PiwsrKyas/j6nOqasmSJdVa\nS2WV6X25/P392bRp09iLFy9YVlYWGzVqlKilj4vPiuvv+WeffcaOHz/OhEIhu3DhAuvTpw9LS0uT\nyfeciluVW0EePHjApkyZwiwtLZmdnR374YcfROvKy8tZaGgos7GxYVZWViw4OJgJhUKp5Ni6dSsz\nNjZmZmZmrHfv3qx3797MzMyMbdy4kRUWFrLZs2eLMh05coTTTN7e3szY2Jj16tVLtMzMzIx5eXlx\nlmnjxo1i271d3GSZ6X25Xr9+zb7++mvRe69du5aVl5fLPNe7MnH1PWeMsaSkJObq6sosLCzYmDFj\n2K+//soYY6ygoEDq33MarJIQopTomhshRClRcSOEKCUqboQQpUTFjRCilKi4EUKUEhU3QohSouJG\nCFFKVNwIAGDp0qUwNjZGt27dYGxsLPbo1q0btmzZwnVEqcjIyMDFixfl/r4//PADduzYAQDw9PTE\n999/L7b+1atXGD16NEaMGIG8vDzcvn0bU6dOlXtOZaLGdQCiGJYtW4YFCxYAAFJSUhAQEICEhATw\n+XwAqDaMUEO1ZMkSWFpawtbWVm7v+eTJExw9ehQnT56scX1JSQmmT5+ON2/eIDo6Gi1atECLFi2g\no6ODEydOYNSoUXLLqkzoyI0AAHR0dESDATZr1gwA0KJFC9GyysEMGzouOuRERETA2dlZ9IeiqrKy\nMvj5+eGff/7B7t27oa+vL1rn7e2NsLAweUZVKlTcSJ0kJyfDw8MDpqamGD58OKKjo0XrYmNj4enp\niR07dsDKygpWVlbYsWMHbt68CScnJ5iZmWHmzJl48+aNaHsPDw9s2bIFn332GWxsbPC///1P7P3O\nnTsHJycn9OrVCy4uLjh16pRo3cKFC7F06VK4urqib9++uHXrFp4/f445c+bAysoKPXv2xPDhw0VD\nDC1cuBApKSnYunUrpkyZgidPnsDY2BgZGRmi14yOjhYN3Hj16lXY2tpi9erV6NOnD9avX//eTG97\n+fIl4uLiqo0ODFQU2oULF+LevXvYvXs3Pv74Y7H15ubmKCgoQGJiokT/NkQcnZYSiWVnZ2P69OkI\nCAjAd999h/v37yMwMBDq6upwd3cHAPz+++9o164dYmNjcerUKXz33Xfo2rUrVq1aBaFQiFmzZiEm\nJgaTJk0CANy9exeamprYv38/Hj58iKVLl6JZs2YYN24c/vjjDyxYsACBgYGwsLBAcnIyli1bhmbN\nmqFfv34AgOPHj2Pjxo0wMDBAly5dMHnyZGhra2Pfvn1QU1PDtm3bsGzZMtja2iIwMBCPHz+Gubk5\nfH19kZ+fDx6PV+33rLrs2bNnyM/Px9GjRwFUDGX+vkxVXb9+Hdra2jA2Nq62LjAwED/99BM2bNiA\nTz75pNp6FRUVWFtb49KlSzW+Nnk3Km5EYlFRUejbty8mTpwIADAwMMCzZ88QGRkpKm5CoRCBgYHQ\n0dHBhAkTsHHjRnh6esLc3BwA0LdvXzx48EDsdb/77jvo6+ujc+fOuHv3Lg4ePIhx48YhIiIC7u7u\nGD16tOj9/vrrL+zcuVO0sxsbG8PR0VH0WkOGDMHQoUPRpk0bAMCUKVNw5MgRPH/+HO3atYOamhq0\ntLTw0UcfIT8//72/M4/Hw4wZM2BgYAAAmD9/fo2ZIiMjayxAt2/fFhsiu1JsbCwEAgG6deuG8PBw\nDB48uMbTViMjIzpy+0BU3IjEHjx4gIsXL4oN5S0UCkUPANDV1YWOjg4AiK7TtW/fXrS9hoaG2OCN\nBgYGYteZTExMsHPnTgDAw4cPcebMGcTGxoq9X+UcBG+/NlDREhkfH4+bN28iPT0dd+7cET3vQ1V9\nDzitDqAAAAMvSURBVEkyVZWXlycapbiqoqIi7Ny5E02bNoWLiws2bdqE+fPnV9uuefPmyM3N/eDs\njRkVNyKx8vJyODk5YdasWdXWqahUXL6tHNG1qppO/SqpqqpWe4/KZeXl5ZgyZQrGjBlT43sBEBtb\nXygUYtKkSSgsLMSIESNgY2MDXV1deHp61vjeNeWqOkJtTe8hSaa336OmwjphwgSYmpoCAPz8/BAa\nGgoHBwf07t1bbDuhUFjra5N3o0+NSKxTp05IT0+HgYGB6PHbb78hIiLig18zMzMTr169Ev2cmpqK\nrl27it4vIyND7P3Onj0ruv71ttu3byMpKQm7du2Cj48P7O3tkZeXB+C/VtKqBY3P54MxJvb+VRsX\nalLXTPr6+jWe/lYt6lOnTkW3bt2wePFilJSUiG2Xn59f64RG5N2ouBGJeXt7488//8S3336LR48e\nISEhAatXr67XzldUVITly5fj4cOHOHHiBA4cOIAvvvgCQMVO//PPP2PHjh148uQJjh07hg0bNojN\n2VCVvr4+VFVV8eOPP+Lvv/8W5QMgOhVu0qQJHj9+jLy8PLRu3RqtWrVCREQEMjIycOrUqVrvRatU\n10w9evTA/fv33/maKioqCA4OxtOnT8Wm4wMqGjDenoKSSIaKG5HYxx9/jG3btiElJQXOzs5YtWoV\nvL294efnV+tz3nVKCgCtW7dGhw4d4O7ujg0bNmDp0qWiiYRNTU2xceNGHDt2DKNGjUJYWBiWLFkC\nV1fXWvOtXLkSe/bswciRIxEaGorZs2dDT08Pt2/fBgB4eHjgl19+wfTp08Hj8bB27Vrcu3cPo0aN\nQkxMTK0zaFWqayZra2uUlJSIFbiaPpMuXbpg5syZOHDgAK5evQqg4pT05s2bGDRo0DszkZrRMOOE\nM7GxsQgLC8OFCxe4jiJTK1asQNOmTUU9QCR18eJFrFu37r1Hk6RmdORGiIx9+eWXiI+Pr3Y97X0O\nHDggmp6Q1B0VN0Jk7JNPPsGYMWOwZ88eiZ9z69YtFBcXw8nJSYbJlBudlhJClBIduRFClBIVN0KI\nUqLiRghRSlTcCCFKiYobIUQpUXEjhCil/wfyw7T+hBgtvAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11b6edf60>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.ticker import LogFormatter, LogLocator\n",
    "from matplotlib.ticker import FormatStrFormatter, MaxNLocator\n",
    "import seaborn as sns\n",
    "sns.set(style='ticks', context='talk')\n",
    "\n",
    "@interact(ilat=[0, len(test_ds.lat)-1], ilon=[0, len(test_ds.lon)-1])\n",
    "def plot_column(ilat, ilon):\n",
    "    \n",
    "    interp_np_column = new_data_np.isel(lat=ilat, lon=ilon)\n",
    "    orig_column = test_data.isel(lat=ilat, lon=ilon)\n",
    "    pres_column = pres_levs.isel(lat=ilat, lon=ilon)\n",
    "\n",
    "    fig = plt.figure(figsize=(4, 5))\n",
    "    ax = fig.add_subplot(111)\n",
    "\n",
    "    ax.plot(orig_column, pres_column/100., color='k', \n",
    "            ms=3, marker='o', lw=1)\n",
    "    ax.plot(interp_np_column, new_pres_levs/100.,\n",
    "            linestyle='none', marker='d', color='r')\n",
    "        \n",
    "    ax.semilogy()\n",
    "    ax.set_ylim(1035, 50)\n",
    "    for p in new_pres_levs/100.:\n",
    "        ax.hlines(p, *ax.get_xlim(), linestyle='dashed', lw=1)\n",
    "    \n",
    "    \n",
    "    ax.set_xlabel(\"{} ({})\".format(test_data.long_name, \n",
    "                                   test_data.units))\n",
    "    ax.yaxis.set_major_formatter(\n",
    "        FormatStrFormatter(\"%d\")\n",
    "    )\n",
    "    ax.yaxis.set_major_locator(MaxNLocator(10))\n",
    "    ax.set_ylabel(\"Pressure (hPa)\")\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Altitude Interpolation\n",
    "\n",
    "We now generalize the interpolation to an arbitrary auxiliary coordinate - in this case, the geopotential height field. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.  2.  3.  4.  5.  6.  7.  8.  9.]\n"
     ]
    }
   ],
   "source": [
    "field = 'Q'\n",
    "test_data = test_ds[field]\n",
    "\n",
    "heights = test_ds.Z3*1e-3\n",
    "new_height_levs = np.arange(1, 10.)\n",
    "print(new_height_levs)\n",
    "\n",
    "# Compute interpolated data\n",
    "new_data_np = ma.interp_by_field(test_data, heights, new_height_levs)\n",
    "\n",
    "# As a safety check, copy attributes into interpolated data\n",
    "nd = ma.copy_attrs(test_data, new_data_np)\n",
    "new_data_np.attrs.update(dict(\n",
    "    long_name='altitude', units='km', positive='up',\n",
    "    standard_name='altitude'\n",
    "))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARsAAAFTCAYAAAAJJqggAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVOX+B/DPsAniMiOSmIyIEgKyiBqKoiJ6XRJ3iH6K\nOYoYCJiFKxZw1X5pCtfLFlheIcEF8CeWRtc0tTK0iBS3q9mgjhQKOqjJzpzfH1wmxhlgYA6z+X2/\nXvMKnnPO8zyH43w6Z2bOdzgMwzAghJAuZqDpCTxPLBYjMTERYrFYZ8fQh32gMbSnf30ZQ+vCprKy\nEklJSaisrNTZMfRhH2gM7elfX8bQurAhhOgnChtCiFpQ2BBC1MJI0xN4Xl1dHQCgtLQUxsbGXTJG\nWVmZ9L9dMUZX909jaNcY+rAPQNNzDvjrOcg2jra99X369GmEhIRoehqEvLBSU1MxadIk1vvVujOb\nwYMHAwCysrJgZWWl4dkQ8uIoKyvDokWLpM9Btmld2BgaGgIArKysYG1treHZEPLiaX4Oso1eICaE\nqAWFDSFELShsCCFqQWFDCFELChtCiFpQ2BBC1ILChhCiFlobNh988IGmp0AIYVGHwqa4uBjjx4+X\n/l5fX48tW7ZgzJgxGDNmDN577z3U19ezMrF9+/ax0g8hRDsoHTa5ubkICgpCQ0ODtC0uLg6//fYb\nvv76a5w4cQK3bt3C3r17WZnY4sWLWemHEKIdlLpdITU1FV999RVCQ0PxySefAAAaGhqQnZ2N3Nxc\n9OzZEwCQmJgoE0aq2LRpEyv9EEK0g1JnNn5+fsjLy4Ozs7O07c6dO5BIJLh48SKmTZuGiRMnYu/e\nvXjppZe6bLKEEN2l1JlN37595doqKytRV1eHM2fO4PDhw3j27BlWrFiBXr16KV0iQiwWy9U7ba7b\nQQjRDJFIJPfaK5fLBY/HU6nfTt/1bWJiAoZhsHr1avTo0QM9evTA0qVLkZmZqXTYZGZmIikpqbNT\nIIR0AYFAINcWHh6OiIgIlfrtdNgMGjQIBgYGMlW9Ghoa0JFaXIGBgfD19ZVpKysrU7izhBD1SE9P\nl6slxeVyVe6302HTs2dPTJ48GfHx8YiLi0NVVRUyMjIwd+5cpfvg8Xhyp2ZdVfKQEKIcPp/fJbWk\nVPpQ37Zt22BlZYXXXnsNc+bMgZeXF5YuXcrW3AgheqRDZzYeHh4oKCiQ/t69e3ds3ryZ9UkRQvSP\n1t6uQAjRLxQ2hBC1oLAhhKgFhQ0hRC20NmyoxAQh+kVrw4ZKTBCiX7Q2bKjEBCH6RWvDhkpMEKJf\ntDZsCCH6hcKGEKIWFDaEELWgsCGEqAWFDSFELShsCCFqQWFDCFELChtCiFpQ2BBC1ILChhCiFhQ2\nhBC1oLAhhKiF1oYN1bMhRL9obdhQPRtC9EuHwqa4uBjjx4+Xa2cYBosXL8ZHH33E2sSong0h+kXp\nsMnNzUVQUBAaGhrklu3ZswdFRUWsTozq2RCiX5QKm9TUVGRmZiI0NFRu2X/+8x8cOXIEU6ZMYX1y\nhBD9oVTY+Pn5IS8vD87OzjLtdXV12LBhA7Zu3Yru3bt3yQQJIfpBqbDp27evwvb4+HhMmDAB7u7u\nrE6KEKJ/OvRd3y0VFBTg/PnzyM3N7fTgYrEYlZWVMm1lZWWd7o8QojqRSIT6+nqZNi6XCx6Pp1K/\nnQ6b/Px8iEQijB07FgBQVVUFQ0NDCIVCpKamKtVHZmYmkpKSOjsFQkgXEAgEcm3h4eGIiIhQqd9O\nh83mzZuxefNm6e8bN24Ej8fDunXrlO4jMDAQvr6+Mm1lZWUKd5YQoh7p6emwsrKSaeNyuSr32+mw\nYQOPx5M7NTM2NtbQbAghAMDn82Ftbc16vx0KGw8PDxQUFChc9uGHH7IyIUKIftLa2xUIIfqFwoYQ\nohYUNoQQtaCwIYSohdaGDdWzIUS/aG3YUD0bQvSL1oYN1bMhRL9obdhQPRtC9IvWhg0hRL9Q2BBC\n1ILChhCiFhQ2hBC1oLAhhKgFhQ0hRC0obAghakFhQwhRCwobQohaUNgQQtSCwoYQohYUNoQQtaCw\nIYSohdaGDRXPIkS/aG3YUPEsQvRLh8KmuLgY48ePl/5+//59hIWFYfTo0fDy8sLWrVvlviO4s6h4\nFiH6Remwyc3NRVBQEBoaGqRta9asQf/+/fH999/j6NGjuHz5MlJSUliZGBXPIkS/KBU2qampyMzM\nRGhoqLStvr4e5ubmCA0NhbGxMSwsLDBr1iz88ssvXTZZQojuUips/Pz8kJeXB2dnZ2mbsbExUlNT\nYWFhIW07ffo0HBwc2J8lIUTnKRU2ffv2bXedrVu3oqSkBCtWrFB5UoQQ/WOkage1tbVYu3Ytfv31\nV2RmZqJPnz5KbysWi1FZWSnTVlZWpuqUCCEqEIlEcm/0cLlc8Hg8lfpVKWweP36M5cuXo0ePHsjO\nzkbPnj07tH1mZiaSkpJUmQIhhGUCgUCuLTw8HBERESr1q1LYhIeHw9LSEomJiTA0NOzw9oGBgfD1\n9ZVpKysrU7izhBD1SE9Ph5WVlUwbl8tVud9Oh80vv/yCwsJCdOvWDaNGjQKHwwEADBs2TOkP5PF4\nPLlTM2Nj485OiRDCAj6fD2tra9b77VDYeHh4oKCgAADg7u6O69evsz4hQoh+0trbFQgh+oXChhCi\nFhQ2hBC1oLAhhKgFhQ0hRC20NmyoeBYh+kVrw4aKZxGiX7Q2bKh4FiH6RWvDhopnEaJftDZsCCH6\nhcKGEKIWFDaEELWgsCGEqAWFDSFELShsCCFqQWFDCFELChtCiFpQ2BBC1EJrw6bl1/wSQnSf1oZN\ndXW1pqdACGERhQ0hRC20Nmxqamo0PQVCCIu0Nmyqqqo0PQVCCIs6FDbFxcUYP3689PcnT54gPDwc\no0aNgo+PD3Jzc1mbGF1GEaJflA6b3NxcBAUFybxL9N5778Hc3BwFBQXYtWsXduzYgeLiYlYmFhAQ\ngMjISFb6IoRonlJhk5qaiszMTISGhkrbqqqqcOrUKaxatQrGxsZwdXXFrFmzkJeXx8rE6uvrkZKS\nwkpfhBDNUyps/Pz8kJeXB2dnZ2nb7du3YWxsjAEDBkjbbG1tIRQKWZmYiYkJVq5cyUpfhBDNU+q7\nvvv27SvXVl1djW7dusm0mZqasvYu0kcffYSwsDBW+iKEaJ5SYaOImZkZ6urqZNpqamrQvXt3pfsQ\ni8WorKyUaSsrKwMA1NbWdnZqhBAViEQi1NfXy7RxuVzweDyV+u102NjY2KC+vh5lZWWwsrICAJSU\nlGDIkCFK95GZmYmkpCSFyyhsCNEMgUAg1xYeHo6IiAiV+u102Jibm8PHxwdxcXHYsmULbt68iWPH\njmH37t1K9xEYGAhfX1+ZtrKyMggEAgobQjQkPT1degLRjMvlqtxvp8MGALZs2YKYmBhMnDgR5ubm\nWL9+PVxdXZXensfjyZ2aGRsbA6AzG0I0hc/nw9ramvV+OxQ2Hh4eKCgokP7eu3dv7Nq1i/VJARQ2\nhOgbrb1dgcKGEP1CYUMIUQutDRu665sQ/aK1YUNnNoToFwobQohaUNgQQtRCa8PGrbQUyMnR9DQI\nISzR2rBZIhQCYWFAebmmp0IIYYHWhk3v+vqmoKE7vwnRC1obNlI5OXQ5RYge0P6wAehyihA9oBth\nQ5dThOg83QgbANk5OYiNjVW4LDY2FhwOR+5B69P6tL7y63fVTdXNOAzDMF06Qgfdu3cPkydPximh\nENbN3+RgaQlcvdr0X0JIl5A+906d6pISE7pxZpOcTEFDiI7T+rC56uQE+PtrehqEEBVpb9hYWKCe\ny0VIY6OmZ0IIYYH2hs3mzTBMTUXJn3/i+vXrmp4NIURF2hs2r70Gg4AALFiwgNXvECeEaIb2hs1/\n+fn5UdgQoge0PmzGjh2LBw8e4ObNm5qeCiFEBVofNoaGhpg/fz4OHz6s6akQQlSg9WEDtLiUopsy\nCdFZKodNUVERFixYgJEjR2LGjBk4duwYG/OSMWHCBFTfvYvGkBC6KZMQHaVS2EgkEoSHhyMkJAQ/\n//wztmzZgg0bNuD3339na34Ami6l9vXqBcNHj+imTEJ0lEph8+TJE4jFYtTX1wMAOBwOjI2NYWho\nyMrkpHJyMFIolPmdLqcI0S0qhQ2Xy8X//M//4N1338WwYcOwePFiREdHo1+/fmzNr/UzGbqcIkSn\ndOi7vp/HMAxMTU2RmJiISZMm4dy5c4iMjISTkxOGDh3a7vZisRiVlZUybWVlZbIrrVypOFSaQyg7\nW5VdIIQ8RyQSSa9WmnG5XPB4PNU6ZlTw1VdfMYGBgTJtkZGRzLZt25TaPiEhgbG3t1f4MDIyYgAw\n2QDDtPbw92cYhmFiYmIYAHKPmJgYhePS+rQ+rS+/fmRkZKvPx4SEBIV9dYRK9WzS09ORn5+PQ4cO\nSds2bNgAS0tLREZGtrt9a2c2AoHgr5oa5eXAsGFyZzdPzczQ884dKj1BCEua69mkp6fDyspKZhkb\nZzYqvWYzduxYXL9+HUeOHAEA/Pjjjzh58iRmzJih1PY8Hg+2trYyDz6fL7uSpWVTPZvnhDQ0UNAQ\n0gX4fL7c81LlSyioGDb29vZISEhARkYGRo0ahS1btmD79u1wcnJSeWIy/P0BPz/przkcDqwiItgd\ngxDSpVR6gRgAvL294e3tzcJU2pGSApw9i9q6Oux3dsaRuLiuH5MQwhqVw0Zt/ns5VSUW47uoKDAM\nAw6Ho+lZEUKUpBP3Rkn5+4MbHAwOh8P6p5QJIV1Lt8IGTZ9SdnNzQ3FxsaanQgjpAJ0LGwBwc3PD\npUuXND0NQkgH6GTYuLq6UtgQomN0MmzoMooQ3aOTYePo6Ijhv/6KuqwsTU+FEKIknQybbk+eIEEi\nAScigu78JkRH6GTYYOVKWDQ2wlgspkJahOgI3QubnByg5Ve7UCEtQnSCboUNFdIiRGfpVti0V0iL\nEKK1tDZsbG1tweFwZB5Xr11rdf3snBy59WNjYxWuGxsbK7curU/rv+jr79q1S+H6bFGpeFZXaC7g\nIy2e1VIrhbRgaQlcvUr1bQhRQZvPPRZo7ZmNQq0U0kJyMgUNIVpOt8IGkCukVcDnN7URQrSa7oUN\n0FRIy9IST7p1wxyRSKl6x4QQzdLNsPnv5VRIYyPKASQmJmp6RoSQduhm2ACAvz/6r1qFbt26oXv3\n7ti0aRMkEommZ0UIaYXuhg2AuLg41NTU4Ndff8Xp06excOFC1NTUaHpahBAFdDpsmllaWuLUqVNg\nGAaTJ09GOX2amBCtoxdhAwBmZmY4cOAAJk6cCE9PT9y4cUPTUyKEtKA3YQMABgYG+N///V9ERUVh\nwoQJOHPmjPxKdOMmIRqhctjcv38fISEhGDlyJLy9vbFv3z425qWSZcuWYf/+/QgICMBnn33214Lm\ne6joxk1C1E7lsFm5ciXs7Ozw008/Yc+ePUhKSsLFixfZmJtKJk+ejNOnTyM2NhbR0dFgGOavGznp\nxk1C1E6lsLl06RLKy8sRGRkJAwMDDBkyBIcOHYKtrS1b81OJk5MTCgoKcOLECSRMmEB1cAjRIJXC\n5urVq7Czs8NHH30ELy8vTJ8+HRcvXkTv3r3Zmp/K+vXrh9PZ2RD89JP8QrqcIkRtVPr63cePH+PC\nhQvw9PTEmTNncPnyZSxfvhx8Ph8jR45sd3uxWIzKykqZtrKyMlWmpJBZZCTMamvlFzRfTmVnsz4m\nIbpKJBKhvr5epo3L5YLH46nWMaOCTz/9lBk7dqxM24YNG5jt27crtX1CQgJjb2+v8GFkZMQAkHnE\nxMQo7CcmJkZuXZn1/fwYBlD4OMRG/7Q+ra8H60dGRrb6fExISFDYV0eoVM/mm2++wcaNG3H+/Hlw\nOBwAwPr169GnTx+sX7++3e1bO7MRCATs1tRopQ5OTc+eMP3tNypPQQj+qmeTnp4OKysrmWVsnNmo\n9JrNuHHjYGZmhqSkJDQ2NqKoqAgnT57EjBkzlNqex+PB1tZW5sHn81WZkmKt1MF5t1s3/D0lBSrk\nLSF6h8/nyz0vVb6Egoph061bN+zbtw+XLl3C2LFjsXbtWrz//vtwdXVVeWKse64ODvz9EX35Mo4f\nP44333wTtYpe0yGEsEalF4iBphT89NNP2ZhL10tJAc6ebfo5ORlWlpY4c+YMAgMDMXXqVBw5cgR9\n+vTR7BwJ0VN6dbtCu5ovp1qUEe3evTtyc3Ph4eEBT09P3Lp1S8OTJEQ/qXxmo3MUlBA1MDDAjh07\nMGTIEHh5eeHw4cMYN26cBiZHiP56sc5s2hESEoKMjAzMmzcPBw4c0PR0CNErL96ZTTumTZuGU6dO\nwdfXF0KhEFFRUdK39QkhnUdnNgq4uLjg/PnzOHLkCJydnWFqakpF1QlREYVNK/r374+zZ8/iP//5\nD2pra5GSkqLpKRGi0yhs2mBubo7ly5cDABYsWNDUSHeLE9IpFDbtSEtLQ2ZmJi5evIjae/eo+BYh\nnURho4SFCxdi6NChuDFlChXfIqSTKGyUwOFw8K8ZM+Dasog6XU4R0iEUNsooL0fvqCj5drqcIkRp\nFDbKaK5d/Dy6nCJEaVobNra2tuBwODKP2NhYhevGxsbKrcvm+levXevQ3Lt6PrQ+rd8V6+/atUvh\n+mxRqXhWV2gu4MNq8SxVtVJ8S9K3LwyuXaPiW0QvdPVzT2vPbLRKK8W30lxcKGgIURKFjbKeK75V\nP28edty+jWPHjmlwUoToDgqbjkhJaTqTsbSEcVoa9u7di7feegsPHz7U9MwI0XoUNh3xXPGtiRMn\nIiAgAGH0jhQh7aKw6Sh/f5kCXB988AEuXryIQ4cOaXBShGg/ChsVmZmZ4bPPPsOqVau65Av2CNEX\nFDYs8PDwQHBwMIKDg+lrYQhpBYUNS6Kjo3H37l1kZGRoeiqEaCUKG5aYmJhg3759WLt2Le7evdvU\nSDdrEiLFWthUVFRg7NixONv8vUwvIFdXV7z77rsICgqC5P59qn1DSAushc2mTZvw+PFjtrrTWWvX\nrsWVK1fwf/37U+0bQlpgJWwOHjwIc3NzuS8jfxEZGRnBu6ICfi1fKKbLKUJUD5uSkhLs3bsXsbGx\n9E4MAJSXY7exsXw7XU6RF5xK3xvV2NiI9evX4/3330evXr06vL1YLEZlZaVMm85/VmXlSvSsrpZv\nb76cys5W/5wI6QCRSIT6+nqZNi6XCx6Pp1K/Kp3ZJCcnw9HREV5eXp3aPjMzE9OnT5d5CAQCANpV\nz6aj67fm6tWrOjF/Wv/FXL+5no1AIJB7XmZmZirsqyNUqmczY8YMVFRUgMNp+sbIp0+fwszMDKGh\noQgODm53+9bObAQCgXbVs+mIVmrfwNISuHqVSlIQrdVczyY9PV3u9Vc2zmxUuozKz8+X+d3Hxwcx\nMTGYOHGiUtvzeDy5HTBW9HqHLmm+WfP112Xb/3vzJiHajs/na3/xrOYznBfec7VvLgwcKHPzJiEv\nIpXObJ536tQpNrvTbSkpwNmzkDAMFj95gpN372LgwIGanhUhGkO3K3SV/15OGaSkYHZQUJcXkyZE\n27F6ZkOe899Lp7fHjIGbmxuio6PB5XI1PClCNIPObNSAz+fD19cXqampmp4KIRpDYaMma9asQUJC\nAmprazU9FUI0gsJGTVxdXeHq6oqsrCxNT4UQjaCwUaO1a9di586dkEgkilegGzaJHqOwUSMfHx+Y\nmpriyy+/lF/YfO8U3bBJ9BSFjRpxOBysXbsWO3bskF+4ciXVvyF6jcJGzfz9/XHx4kUYGRkhJCSk\nqTEnB8jN/Wslupwieog+Z6NmRkZGqKmpQWNjI9LS0mBcWYn4Eycgd0dYWBjg7U33UxG9QWc2GhAe\nHg5TU1OEhoYi7OpVGIvF8ivR5RTRM1obNrpcz6a99ePi4lBdXY2UlBQ4ODi0+jeg+je0vibq2XQV\nlerZdIXmmho6W8+mo1qpf1PO4SBu6VIsW78e9vb2GpoceZF09XNPa89sXhjN9W+eY/jxxzDq3x9e\nXl6YOXMmTpw4QTWeiU6jsNEGz9W/gb8/+rz1FrZu3Yo7d+5g/vz5WLNmDZydnZGWloaqqirNzZWQ\nTqKw0RYpKU1nOc+d6ZiZmSEoKAiXLl1CYmIivvzyS9jY2GDjxo24d++eBidMSMdQ2GiL5pBppXwo\nh8OBj48Pjh49ivPnz6O6uhpubm544403UFBQQJdYROtR2GgTf3+lyocOGTIEu3btQklJCTw9PbF4\n8WKMGTMG+/fvR11dnRomSkjHUdjosF69euHtt9/GjRs3sGnTJuzZswe2trb44IMPUE73VxEtQ2Gj\nBwwNDTF79mycOnUK+fn5EAqFsLe3x/Lly3H58mVNT48QABQ2esfV1RV79uzBzZs3YWtri+nTp2Py\n5Mn4/PPP0djYqOnpkRcYhY2esrS0xKZNm1BSUoKgoCBs3boVQ4cOxaRJk2BqaorIyEhNT5G8YFQO\nm8LCQrz++usYNWoUpk6dikOHDrExL8ISExMTLFy4EBcuXMC+ffvw7bffora2Frt27UJWVhaePHnS\n8U7prnTSCSqFzZMnTxAWFgaBQIDCwkLs2rUL8fHxKCgoYGt+hCUcDgeenp5YvXo1TE1NMXXqVBw4\ncAB8Ph/z5s3D/v378fTp0/Y7oiJfpJNUCpvff/8d3t7eeO211wAATk5OGD16NH755RdWJkfY13wT\naH5+Po4dO4Y7d+5g7ty5yMrKgrW1NebNm4cDBw60HjxU5It0kkph4+DggO3bt0t/f/z4MQoLC+Ho\n6KjyxIh6cLlcLFmyBMePH8ft27cxZ84cZGZmwtraGvPnz8fBgwfx559/Nq1MRb6IClh7gfjp06cI\nCQmBi4sLJk2axFa3RI14PB4EAgGOHz+OkpISzJo1C5999hkGDBiApb6+qFm+XH4jupwiSmKlUp9I\nJEJoaChsbGzwj3/8Q+ntxGIxKisrZdrKysrYmBJRUZ8+fbB06VIsXboUjx49wuOpU2Gq6MXk5sup\n7Gz1T5J0CZFIhPr6epk2LpcLHo+nWseMiq5cucJ4enoy27Zt6/C2CQkJjL29vcKHkZERA0DmERMT\no7CfmJgYuXVpfXbXzwYYppXHOWtrrZ8/rd/++pGRka0+HxMSEhT21REqFc+qqKjA7NmzsWzZMixX\ndIrdjtbObAQCwYtTPEtXtFLkq7pnTywZORL/LirCtGnT4O/vj9deew3m5uYamijprObiWenp6bCy\nspJZxsaZjUqXUYcPH4ZYLEZKSgqS/1sWgcPh4M0338Tq1avb3Z7H48ntgLGxXOlvog2a70p//XWZ\nZrM9e5Dt74+Kigrk5eXhk08+wfLlyzF9+nRp8HTv3l1Dkyadwefzu+Z/9CqfG7FMJBIx9vb2jEgk\n0vRUiCJ+fn9dQvn7K1ylvLyc2b17NzNlyhSmV69eTEBAAJObm8s8e/ZMzZMlHdHVzz36KhfSMSkp\nwNmzTT8rKGcKAH379kVwcDCCg4NRXl6OI0eO4N1334Wfnx94PB5Gjx4NW1tbDB48WPqwtbVF7969\n1bgjRN0obEjHtKwkqMR3WllaWmLFihVYtWoVAKCqqgrh4eEQCoUoKSnB999/D6FQCKFQiG7duskE\nUMsg4vP5dImt4yhsSMcpUeDreWFhYUhJScHKlSsxc+ZMueUMw6C8vFwaQkKhEOfPn8f+/fshFApR\nVlaGAQMGyARQy0Dq06cPOBwOG3tHugh9lQvRCXV1dbh79670LKhlKP32229gGEYugJpDadCgQejW\nrZumd0HrdfVzj85siE4wMTGBnZ0d7OzsFC4Xi8UyIVRcXIy8vDwIhUKIRCK89NJLciHU/HO/fv3o\nrEgNKGyIXuDxeBg5ciRGjhwpt6yhoQGlpaUyZ0XHjx+XBtOzZ89ga2vb6pkRvXXPDgoboveMjIxg\nY2MDGxsbhfftPX36VHpJJhQKcevWLZw4cQJCoRC3b99G7969peFTUlKCn3/+GWFhYYiLi5PtqPmm\n1E68pvUioLAhL7yePXvC1dUVrq6ucsskEgn++OMPaRgtXboUEokEKSkpsmHTsuSGt7dS79S9aKgs\nKCFtMDAwwIABA+Dl5YU333wTwcHBAAC/lt9gClCdHyVQ2BDSAampqfjXv/6Fa9euoaGhoamR6vwo\nhcKGkA4SCASoqKiAmZkZokNDFZ/JUJ0fORQ2hHQQh8NBWVkZGhoa4JKaqjhU6HJKjtaGja2tLTgc\njswjNjZW4bqxsbFy69L6tH5Xrq/M1xxfuXoVeXl5WL16tdbNX9H6u3btanefVEGfICZEFa3U+anp\n1Qvxy5bh3M2bKCoqQl1dHdzd3eHu7o4RI0bA3d0dr7zyCgwNDTU0cXn0CWJCtFkrdX5MP/0UUS0+\nb1NWVoZffvkFRUVFOHz4MDZt2oQHDx7Azc1NJoCGDRsGExMTde+FWlDYEKIqf3/Az++vd6T8/eU+\n2GdlZYUZM2ZgxowZ0rbKykpcvHgRRUVFOH36NOLj43Hz5k1IJBIEBwfj448/VudedDkKG0LYoESd\nn+dxuVx4e3vD29tb2mZqaor6+nqkpqaiR48eWLt2LV566aUumLD6ae0LxITolObLqeRklT49HBYW\nBlNTU6xYsQI1NTVwdHTEunXr8ODBAxYnqxkUNoSwRcHlU0c1f2NpWloaEhMTcenSJVRVVcHR0RGv\nvvoqTE1NERkZydKE1YvChhAtZm1tjaSkJFy8eBE///wzamtrkZKSoulpdQqFDSE6gM/nY8qUKTA0\nNMTKlSs1PZ1OobAhREd88sknMDExQXJysk5eSlHYEKIjbGxsUFNTo7OXUiqHzbVr1+Dv7w93d3fM\nmzcPly5dYmNehBAFXFxcYGRkpJOXUiqFTV1dHUJDQ+Hn54fCwkIEBgYiNDQU1dXVbM2PENLCwIED\nIZFIND2NTlEpbM6fPw9DQ0MEBATA0NAQCxYsgIWFBc42f7iJEMKqr7/+GhKJRPp117pEpbARCoUY\nMmSITJutrS2EQqFKkyKEKBYWFgZDQ0M4OTlpeiodplLYVFdXw8zMTKbNzMwMNTU1Kk2KEKJYXFwc\n7t27h5ty35njAAAVAUlEQVQ3b+rcB/xUujdKUbBUV1cr/dUXYrEYlZWVMm2///47gKa7ZAkhijEM\ng8bGRqSlpeGdd95hpc/m59ydO3dQX18vs4zL5YLH46nUv0phM3jwYGRlZcm0lZSUYPbs2Uptn5mZ\niaSkJIXLFi1apMrUCNFrLevNTJ48mdW+ly1bJtcWHh6OiIgIlfpVKWzGjBmDuro6ZGVlISAgAHl5\neXj06BG8vLyU2j4wMBC+vr4ybUKhECtXrsSePXswaNAgVabXKpFIBIFAgPT0dPD5fJ3rn8bQrjH0\nYR8A4Pbt2wgKCkJKSgoGDx4ss4zL5arcv0phY2Jigk8++QTR0dGIj4+HjY0NPv74Y5iamiq1PY/H\na/XUbMCAAV1Wqa/5FNHKyqpLxujq/mkM7RpDH/ah5RjN3wTKNpXr2djb2+PgwYNszIUQosfodgVC\niFpQ2BBC1MIwtrXve9AgU1NTeHh4yH2GR5fG0Id9oDG0p399GEPrvsqFEKKf6DKKEKIWFDaEELWg\nsCGEqAWFDSFELShsCCFqQWFDCFELChtCiFp0WdgoWwj92LFjmDJlCtzd3RESEoKHDx8q1ceTJ0/w\n5ptvwtHREQ4ODvD29mZ9jPv372Px4sVwcnLC0KFD4enpiZ9//pnVMVouc3d3x5o1a1jtv76+Hu+8\n8450H0aPHs36Pjx48ACLFi2Ck5MTHBwcMG7cuE4di2bFxcUYP368TBtbx7utMdg63m2Nwcbxbqt/\nto53W2M8ePAAISEh8PDwwPjx4/GPf/xDYf9ymC5QW1vLTJgwgTl48CDT0NDA5ObmMp6enkxVVZXM\netevX2dGjhzJFBcXM7W1tcymTZuY4OBgpfoICwtj3NzcmKysLKaoqIhxc3NjXn31VVbHWLhwITN8\n+HAmKyuLuX//PjN58mTGzc2NtTFaLktLS2McHBwYV1dXVvdh69atjLOzM5Oens48evSImTJlCjN8\n+HDWj8Xw4cOZAwcOMKWlpczo0aOZESNGdGiMZjk5OcyoUaOYMWPGyLSzcbzbG4ON493WGGwc7/b2\ngY3j3d4YERERzIcffshIJBKmrKyMmTx5MpOXl8e0p0vObJQthN6crC4uLjAxMcGaNWvw3Xff4dGj\nRygoKGi1j6qqKpw6dQq9evXCwoUL4e7uDj8/P3A4HNbGqK+vR11dnXSMl156CQKBgNUxmv9Obm5u\nOHr0KKZOnQozMzPW+m9oaMChQ4fA5XKxZMkS8Hg87N+/H1ZWVqyNATT935rD4WDBggUwMDBA7969\n0atXrw6NAQCpqanIzMxEaGiozHZsHe+2xmDreLc1BhvHu63+2TrebY0BNBXIa2hoQENDAxiGgaGh\noVJlZbokbJQthP78elwuF1wuF0KhECUlJa32cefOHRgaGmLo0KEyywwMDFgbw9jYGDNnzoS9vb10\n2enTp/HSSy+xNoZQKIStrS02bNiArVu3onv37ujduzerf6fGxkZwuVxMmzYNEydOxN69e2FnZ8fa\nGAAwYsQIVFdXw93dHZMmTcKIESPg4uKi1Bgt99fPzw95eXlwdnaW2U7V463MGKoeb2XGUOV4K/t3\nkkgknT7eyowBAMuXL0d2drbM8Z42bZrces/rkrBRthC6ovVMTU1RU1PTZh9VVVUwNDSUWd6crGyN\n8fy2W7duRUlJCRwdHVkbo7q6GiKRCBMmTIC7uzsAwMjIiLX+Kysr0djYiMePH+Pw4cPIzs7GuXPn\nIBKJWP071dXVYfDgwSgqKsKxY8dQWFiI+/fvKzVGy3769u0LRVQ93sqM8fy2HT3eyoyhyvFWpv/K\nyko0NDR0+ngr+3diGAYhISEyxzs7O7vV9Zt1SdgoWwi9+R+zovXa6sPMzAyNjY0yy2tqasDhcFgb\no3k/qqqqsGrVKpw7dw6ZmZmQSCSsjXH//n1UVFRg1apV0mUNDQ2s9W9iYgKGYTBw4ED06NED/fr1\nw9KlS1FaWsraGOXl5fjmm2/Qv39/mJiYYMiQIVixYgWEQmGHxmgLW8e7PWwc77awcbzbYmJiAgAq\nH++2lJeXIzY2FsHBwTLH+9ChQ21uB3RR2AwePBglJSUybSUlJbCzs5NpGzJkiMx6jx49wpMnTzBk\nyJA2+7CxsYFEIsGtW7dkljU2NrI2BtBUgvHChQt4+vQpsrOz8fLLL7O6H6WlpaiqqsLYsWPh4eGB\nL774Anfu3MHXX3/NSv+DBg2CgYEBRCKRdFlDQwNqa2tZ24fy8nJIJBLcvn1buszAwAB1dXUdGqMt\nbB3v9rBxvNvCxvFuy6BBg8DhcFQ+3m0pLy+Xvl7TzMDAAEZG7Rf97JKwaVkIvaGhAbm5uQoLofv6\n+uLEiRMoKipCbW0t4uPjMWHCBPTu3bvNPszNzeHj44OHDx8iPT0dv/zyCw4fPgyGYVgbAwD27dsH\nDocDHx8fmJmZsb4fycnJ6NevH95++2388MMPcHFxgampqdw3VnS2/549e2Ly5MkQi8XYs2cPSktL\nkZCQACMjI9b2wc7ODlZWVqioqEBGRgbu3LmDXbt2wdDQsENjtIWt490eNo53W9g43m1h63i3xc7O\nDv369cO2bdtQV1eHe/fuYe/evZg5c2ab2wHomre+GYZhbty4wQQEBDAjRoxg5s2bx1y6dIlhGIaJ\njo5mYmJipOvl5+czU6dOZUaOHMm89dZbzMOHD9vtg2EYprKyklm2bBnj6OjIDB06lJk0aRKrYxQV\nFTEODg6Mi4sL4+DgwNjb2zMODg7M3LlzWd2Plss8PDyYNWvWsNr/s2fPmIiICGbYsGGMvb094+np\nyVy8eJHVMW7dusUEBARIj4WXl1enjkWzCxcuyL3dytbxbm0MNo93W/vBxvFuq3+2jndbY9y6dYtZ\ntmwZ4+HhwUyaNIn5+OOP5bZThIpnEULUgm5XIISoBYUNIUQtKGwIIWpBYUMIUQsKG0KIWlDYEELU\ngsKGEKIWFDbtaGxsRGpqKqZNmwYXFxeMGzcO69evx++//66W8ZOSkvDGG29If8/IyICHhwdeffVV\nnDx5Eo6Ojqirq+twvz4+Pkrdz8IWBwcHfP/99wqXlZaWwsHBQfrx+cWLFyM+Ph5AU+mHAwcOqDT2\n7t27sXv3brm+OysiIgInT55EUlISAgICOtVHeXk55s+fL/Oxf31HYdOOuLg4HD16FNHR0fj3v/+N\njz/+GBUVFQgMDERtbW2Xjx8UFCR9ogBAfHw8BAIBPv/8c3h7e+P777+X3oCnzc6dO4cxY8a0upzD\n4Uh/Tk5OltZROX78OJKTkzs97r1795CTkwOBQNDpPlqSSCQoLCyU7kvLeXeEpaUlvLy8kJaWxsq8\ndAGFTTsOHz6MiIgIjBs3Di+//DJcXV3xz3/+E+Xl5XIFibqCmZkZevXqBaCplENtbS1GjRqF/v37\nw8jICBYWFl0+BzZYWFi0ebNeyw+y9+rVS1r+QCKRqDTunj174Ovry1ogX758GTY2NujRo4fKfS1a\ntAgZGRmorq5mYWbaj8KmHQYGBigoKEBjY6O0rUePHjh27Jj05rakpCRERERg8+bNcHd3h4+Pj1x9\nj8OHD2PatGkYPnw4/P39UVBQIF3GMAySk5Ph7e0Nd3d3BAUF4e7du9K+AwICUFpaCldXV3A4HCxZ\nsgQbN27Ejz/+CAcHB+ll1B9//IGVK1dixIgR8PLyws6dO9HW3SglJSVYvHgxXF1dMWPGDJw7d066\n7PnLnu+++w4ODg4yy/Pz8zF79my4ubkhKCgIZWVlWLduHdzd3TFt2jRcuHBBYX/V1dXYuHEjRo0a\nBR8fH3z77bcyZwjNlzo//vgjoqKiUFFRAUdHRxQXF8PBwUHmEvbx48dwdnbG9evX5fbv2bNnyMvL\nw5QpUxTuf0VFBaZOnYp169ZJ2zIzMzFp0iS4u7tj3bp1iIyMRFJSknT5uXPn4OnpKdeXRCLBypUr\nMXv2bDx58gQAUFBQgLlz58LNzQ2BgYFITEzE4sWLpdv069cPAwcORF5ensL56RsKm3YsXboUOTk5\nmDhxIqKionD06FE8evQINjY2MrU/zpw5A7FYjNzcXISGhmLz5s04c+YMAODs2bPYuXMn1q5diy++\n+AJz585FSEgIfv31VwBAYmIi9u/fj5iYGBw9ehTm5uYICwuT9s3hcPDyyy/jm2++AcMwSExMxKZN\nm6TLgKazHoFAgMbGRhw6dAi7du1CXl4e9u7d2+q+5eTkYNGiRTh+/DgcHR1lnnSKPH/JEB8fj5iY\nGGRmZuLKlSuYPXs2nJyc8H//93+ws7PD5s2bFfbz/vvv48qVK9i7dy927tyJjIwMheuNGDECUVFR\nsLCwwLlz5+Ds7AwbGxvk5+dL1/n3v/8NGxsbODo6ym3/448/olu3bhg2bJjcsj///BPBwcFwcHDA\ntm3bADRdsu3cuRPvvvsuDh8+DENDQ3z55Zcy2xUUFGDs2LFy/b333nsQCoVIT09Hr169cO/ePYSE\nhMDb2xtHjx7F3/72N6Slpcn9Db28vPDtt98q3H99Q2HTjhUrVuCf//wn7Ozs8Pnnn2PDhg2YMGEC\ntm/fLrOeubk5tm3bhiFDhsDf3x8zZ87EwYMHAQCffPIJgoODMWXKFPD5fCxatAhTpkzBZ599BgA4\nePAgIiIiMGnSJAwcOBAxMTGYMGECqqqqpP1zOBxp9bTevXvLncb/8MMPuH//Pnbs2IFXXnkFo0aN\nwubNm9GnT59W9+3111/H9OnTwefzsWLFCjx69AhlZWVK/20WL16MkSNHwsXFBZ6enhg4cCAEAgFs\nbW0REBCA27dvy51Z/fnnn8jPz0dUVBRcXFwwYsQIbNy4UWH/RkZG6NmzJzgcDvr06QMDAwPMmjUL\nx48fl65z/Phx+Pr6Ktz+6tWrcnVcgKbCW6GhobC0tER8fDwMDJqeBvv27cOiRYswa9YsDB48GFu2\nbEG/fv2k21VVVeHmzZvSKnvNdu7cifPnzyMjI0P6987JycErr7yC1atXY9CgQViyZAmmTp0qNxc7\nOztcvnxZ4fz1TfsVbwimTZuGadOmoaqqCufPn0deXh7S09MxYMAABAYGAgCcnJzQrVs36TYuLi7S\nMLl16xYuX76MxMRE6fKGhgYMHz4cjx49wqNHj2RqvVpYWGDt2rUdmuNvv/0GPp8vfX0HaHrHqS18\nPl/6c8+ePcEwjFz1trYMHDhQ+rOZmZlMf6amppBIJKivr5d5vaSkpAQSiUTmkszFxaXNy72WZs2a\nheTkZNy9exempqYoLCzE1q1bFa778OFD8Hg8ufaDBw+ioaEBc+fOlXkd6caNG1i6dKn0dyMjI5nj\n8tNPP8HNzU1mm2vXrqG4uBiDBg2Sef3sxo0bcHNzkxl3+PDhOHnypEwbl8tFZWWlUvuu6+jMpg03\nbtyQ+YfcvXt3+Pj4ICEhAT4+PjKvaRgaGsps29jYKP0/ZmNjI9atW4fPP/9c+vjyyy+xc+dO1l64\nNDY27vA2z88ZQKtP+pavWTV7/gXf5v1tS/NlRMtxOjJ3GxsbuLq6Ij8/HydOnICLi4tMyD0/H0Uv\nMNvZ2SEtLQ1HjhxBYWGhzDzaekH63LlzcpdQpqamyMzMxOPHj2XeNTQyMlLqxe3GxsZOv6Olayhs\n2iCRSJCZmSnzD7JZjx49ZC5Rbt68KfOPq/nFTKCpBGNpaSn4fL70kZOTg6+//ho9evSAhYWFzAuc\njx8/hqenp1xF/LYMGjQI9+7dw59//ilty8rKQlBQUIf2uZmxsbFMX80vWKvK1tYWRkZGKC4ulrZd\nuXKl1SeconZfX1988803OH36NGbNmtXqWH379oVYLJZr9/Lywvjx4zFnzhzExsZKg/SVV17B1atX\npetJJBJcu3ZN+rui12sGDx6MESNGIDIyEmlpadKSnM/3BUDh5ZJYLNaZdxRVRWHTBkdHR0ydOhWr\nV69Gbm4uRCIRrl+/jt27d+PkyZNYsmSJdN0HDx5IK/Lv378fJ06ckL7zsHz5cuzfvx85OTkQiURI\nT0/Hnj17YGtrCwBYsmQJkpKS8N1330EoFOL9999H//79MXjw4Hbn2HyG4OXlBWtra0RFReHWrVso\nKChAWloaJkyY0Kl9d3FxQVZWFkpKSvDDDz9g3759He5D0VmSubk5FixYgA8//BCFhYW4dOmS9AVa\nRbp3745nz55BKBRKQ+G1117D1atXUVRUhBkzZrS67bBhw6Qvwiuydu1alJWV4dNPPwUACAQC7N+/\nH8eOHUNJSQk2b96MP/74AxwOBw8ePIBYLJb5qpeW5s+fD0dHR8TGxgIA3njjDdy8eRMJCQm4ffs2\nDhw4gPz8fLnwvHHjhsKvS9FHFDbtiIuLw8KFC5GRkYHZs2cjMDAQP/74I9LT02W+x8jR0RGNjY2Y\nN28e9u3bh/j4eIwcORIAMGXKFGzatAmffvopZs6ciZycHOzYsQPjxo0D0PTBvblz52Ljxo3w9/dH\nXV1dqx9ke/4fa/PvBgYGSElJQU1NDfz9/bFhwwa8/vrrMoHYVj/Pt0VHR6O2thZz5sxBXFwc3nnn\nnXa3b6u/lj9HRUVh/PjxCA0NRXh4OBYuXNjqdmPGjMGQIUMwb9486VmGhYUFRo8ejVGjRrX5AviY\nMWPQ0NAgc9bYsm8LCwtEREQgNTUVpaWl+Nvf/oaIiAhs374d8+fPR319PYYPHw5jY2P88MMPbX4o\nEWh6l+38+fP48ssv0a9fP6SkpODEiROYPXs28vPzMWfOHLlLxsLCQkycOLHNfvWGUsVDSZsSExOZ\ngIAATU/jhTJr1izmiy++aHe9v//978y2bduU6vPChQvMvXv3ZNqmT5+u1FfLPu/mzZvMlStXZNqi\no6OZqKgo6e8lJSWMh4cH8+zZsw73r4vozIbolO+++w7x8fF4+PChwreSn7ds2TJ89dVXSn1K95tv\nvkFYWBiuXLkCkUiE5ORkPHz4EOPHj+/wPEUiEd58802cPXsWv//+O/Lz83Hs2DGZbyE4ePAgFi9e\nrNR3WukDeuub6JR9+/bhypUr+PDDD5V6J8/a2hpvvPEGMjIyEBIS0ua6b7/9Np4+fYoVK1aguroa\nTk5O2LNnT5uXaq3x8fFBSEgI/v73v+Phw4ewtrZGdHS09AXmBw8e4Pz580p9k6S+oG9XIISoBV1G\nEULUgsKGEKIWFDaEELWgsCGEqAWFDSFELShsCCFq8f8E855bJwejwQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11c47d7f0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@interact(ilat=[0, len(test_ds.lat)-1], ilon=[0, len(test_ds.lon)-1])\n",
    "def plot_column(ilat, ilon):\n",
    "    \n",
    "    interp_np_column = new_data_np.isel(lat=ilat, lon=ilon)\n",
    "    orig_column = test_data.isel(lat=ilat, lon=ilon)\n",
    "    heights_column = heights.isel(lat=ilat, lon=ilon)\n",
    "\n",
    "    fig = plt.figure(figsize=(4, 5))\n",
    "    ax = fig.add_subplot(111)\n",
    "\n",
    "    ax.plot(orig_column, heights_column, color='k', \n",
    "            ms=3, marker='o', lw=1)\n",
    "    ax.plot(interp_np_column, new_height_levs,\n",
    "            linestyle='none', marker='d', color='r')\n",
    "\n",
    "    for p in new_height_levs:\n",
    "        ax.hlines(p, *ax.get_xlim(), linestyle='dashed', lw=1)\n",
    "    \n",
    "    ax.set_ylim(0, 17)\n",
    "    ax.set_xlabel(\"{} ({})\".format(test_data.long_name, \n",
    "                                   test_data.units))\n",
    "    ax.yaxis.set_major_formatter(\n",
    "        FormatStrFormatter(\"%d\")\n",
    "    )\n",
    "#     ax.set_ylabel(\"{} ({})\".format(new_data_np.lev.long_name, \n",
    "#                                    new_data_np.lev.units))\n",
    "    "
   ]
  }
 ],
 "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.4.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}