Read_Grid_Interpolate.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import pygrib\n",
    "import numpy as np\n",
    "from netCDF4 import Dataset\n",
    "import os\n",
    "from scipy import spatial\n",
    "import pandas as pd\n",
    "def lon_convert(station_lon):\n",
    "    if station_lon<0:\n",
    "        diff_from_180=180+station_lon\n",
    "        station_lon_360=180+diff_from_180\n",
    "        station_lon_180=station_lon\n",
    "    else:\n",
    "        diff_from_360=360-station_lon\n",
    "        station_lon_180=-diff_from_360\n",
    "        station_lon_360=station_lon\n",
    "    return(station_lon_180,station_lon_360)\n",
    "\n",
    "#This came from Stack Overflow! Thanks random internet person\n",
    "#https://stackoverflow.com/questions/8661537/how-to-perform-bilinear-interpolation-in-python\n",
    "def bilinear_interpolation(x,y,coord_vals):\n",
    "    #We are assuming that the coord_vals are a list of 4 tuples containing\n",
    "    #x,y,and value\n",
    "    coords=sorted(coord_vals)\n",
    "    (x1,y1,q11),(_x1,y2,q12),(x2,_y1,q21),(_x2,_y2,q22) = coords\n",
    "    if x1 != _x1 or x2 != _x2 or y1 != _y1 or y2 != _y2:\n",
    "        raise ValueError('points do not form a rectangle')\n",
    "    if not x1 <= x <= x2 or not y1 <= y <= y2:\n",
    "        raise ValueError('(x, y) not within the rectangle')\n",
    "        \n",
    "    return (q11 * (x2 - x) * (y2 - y) +\n",
    "        q21 * (x - x1) * (y2 - y) +\n",
    "        q12 * (x2 - x) * (y - y1) +\n",
    "        q22 * (x - x1) * (y - y1)\n",
    "       ) / ((x2 - x1) * (y2 - y1) + 0.0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def interpolate_around_station(ystart,yend,var_string,lat_south,\n",
    "                              lat_north,lon_west,lon_east,\n",
    "                               station_lat,station_lon):\n",
    "    #Convert longitude coordinates\n",
    "    station_lon_180,station_lon_360=lon_convert(station_lon)\n",
    "    mon_str=[\"{num:02d}\".format(num=m) for m in range(1,13)]\n",
    "\n",
    "    print('reading in variable '+var_string)\n",
    "    #Find the points (read in lat/lon data from example grib file)\n",
    "    fstr='/gpfs/fs1/collections/rda/data/ds627.1/ei.moda.an.sfc/ei.moda.an.sfc.regn128sc.1979010100'\n",
    "    grib_open=pygrib.open(fstr)\n",
    "    grib_var=grib_open.select(name=var_string)[0]\n",
    "    data,lats,lons=grib_var.data(lat1=lat_south,lat2=lat_north,\n",
    "                                       lon1=lon_west,lon2=lon_east)\n",
    "    #make the lat and lon axes\n",
    "    lat_axis=lats[:,0]\n",
    "    lon_axis=lons[0,:]\n",
    "        \n",
    "    pt=[station_lon_360,station_lat]\n",
    "\n",
    "    #Find the 4 nearest coordinates\n",
    "    #First in the x direction\n",
    "    x_dist=[abs(station_lon_360-l) for l in lon_axis]\n",
    "    #Then the y\n",
    "    y_dist=[abs(station_lat-l) for l in lat_axis]\n",
    "\n",
    "    idx=np.argpartition(x_dist,2)\n",
    "    idy=np.argpartition(y_dist,2)\n",
    "\n",
    "    nearest_inds=[]\n",
    "    nearest_pts=[]\n",
    "\n",
    "    #make tuples of the coordinate pairs\n",
    "    for i in range(0,2):\n",
    "        for j in range(0,2):\n",
    "            nearest_inds.append((idx[i],idy[j]))\n",
    "            nearest_pts.append((lon_axis[idx[i]],lat_axis[idy[j]]))\n",
    "            \n",
    "    interp_vals=[]\n",
    "    dates=[]\n",
    "    mon_sum=np.zeros(12,dtype=float)\n",
    "    for y in range(ystart,(yend+1)):\n",
    "        mstart=0\n",
    "        mend=12\n",
    "        if y==ystart:\n",
    "            mstart=mon_str.index('10')\n",
    "        elif y==(yend):\n",
    "            mend=mon_str.index('09')+1\n",
    "\n",
    "        for m in range(mstart,mend):\n",
    "            fstr='/gpfs/fs1/collections/rda/data/ds627.1/ei.moda.an.sfc/ei.moda.an.sfc.regn128sc.'+str(y)+mon_str[m]+'0100'\n",
    "            #print('opening '+ fstr)\n",
    "            #Open the grib file\n",
    "            strdate=str(y)+'-'+mon_str[m]\n",
    "            dates.append(strdate)\n",
    "            #print(strdate)\n",
    "            grib_open=pygrib.open(fstr)\n",
    "            grib_var=grib_open.select(name=var_string)[0]\n",
    "            #subset the variable\n",
    "            var,lats,lons=grib_var.data(lat1=lat_south,lat2=lat_north,\n",
    "                                       lon1=lon_west,lon2=lon_east)\n",
    "\n",
    "            #grab the value that corresponds to the lat/lon pairs\n",
    "            nearest_copy=nearest_pts.copy()\n",
    "            for i in range(0,len(nearest_inds)):\n",
    "                val=var[nearest_inds[i]]\n",
    "                nearest_copy[i]=nearest_copy[i]+(val,)\n",
    "\n",
    "            new_val=bilinear_interpolation(station_lon_360,\n",
    "                                           station_lat,nearest_copy)\n",
    "            mon_sum[m]+=new_val\n",
    "            interp_vals.append(new_val)\n",
    "            #close the grib file\n",
    "            grib_open.close()\n",
    "    ny=yend-ystart\n",
    "    mon_avg=mon_sum/ny\n",
    "    \n",
    "    interp_copy=interp_vals.copy()\n",
    "    anom_vals=[]\n",
    "    for i in range(0,len(interp_copy)):\n",
    "        date_month=dates[i][5:]\n",
    "        mon_ind=int(date_month)-1\n",
    "        anom_vals.append(interp_copy[i]-mon_avg[mon_ind])\n",
    "    \n",
    "    return(interp_vals,anom_vals,dates)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "reading in variable 2 metre temperature\n"
     ]
    }
   ],
   "source": [
    "ystart=1979\n",
    "yend=2017\n",
    "var_string='2 metre temperature'\n",
    "lat_south=32\n",
    "lat_north=52\n",
    "lon_west=234\n",
    "lon_east=260\n",
    "station_lat=38.5050\n",
    "station_lon=-107.8975\n",
    "\n",
    "temps,anoms,dates=interpolate_around_station(ystart,yend,var_string,lat_south,lat_north,lon_west,lon_east,station_lat,station_lon)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "index = pd.date_range(start='10-01-'+str(ystart), end='09-30-'+str(yend), freq='M')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame(data=anoms, index=index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x2ae2ca568780>"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEHCAYAAABMRSrcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd83GeV6P/PmVEbFat3ySW2LHcnsWOn2I5DDGkkWSCEFC6wAUIgsLuU3cte9l728gN29wJ7CYRLEloWNiSEEmKTBgk4dpptJRv3GjcVW71aozbz/P6YGXniSLKkKd/vd+a8Xy+/LI1m5vs4+UpHz3POcx4xxqCUUiq5uawegFJKKetpMFBKKaXBQCmllAYDpZRSaDBQSimFBgOllFJoMFBKKYUGA6WUUmgwUEopBaRYPYCJFBUVmdmzZ1s9DKWUcozXX3+9zRhTPNXX2ToYzJ49m7q6OquHoZRSjiEiJ6bzOl0mUkoppcFAKaWUBgOllFLYPGeglFJWGB4epqGhgYGBAauHMq6MjAyqqqpITU2NyvtpMFBKqXM0NDSQk5PD7NmzERGrh/MOxhja29tpaGhgzpw5UXlPXSZSSqlzDAwMUFhYaMtAACAiFBYWRnXmosFAKWULIz4/djp50a6BICTa49NgoJSyhc888gbX3beVlh77rtMnMg0GSinL9QwM88KBFg6c7uXWB1+lsctr9ZBs4dlnn6W2tpZ58+bxr//6rzG9lgYDpZTlXj7chs9v+MfrFtB+ZohbH3iV421nrB6WpXw+H/feey/PPPMM+/bt49FHH2Xfvn0xu55WEymlLLf5YCs5GSnctWYOV8wr4iM/3c4HH3yVRz6xmvmlOZaO7X9v2su+pp6ovueiihl89cbFEz5n+/btzJs3jwsuuACA2267jSeffJJFixZFdSwhOjNQSlnKGMOLh1pZW1NEqtvFkspcfnX3pQjwoQdfZU9jt9VDtERjYyPV1dWjn1dVVdHY2Biz60VlZiAiPwXeC7QYY5aM8fX1wJPAseBDvzPGfC0a11ZKOduB072c7hlg/fyS0cdqSnP49T2XccePtnH7j17j4b9exYpZ+ZaM73y/wSeKaM0MHgauPc9zthpjLgz+0UCglAICS0QAV9a+vevyrMIsHr/nMoqy0/lvP9nGK2+1WTE8y1RWVlJfXz/6eUNDA5WVlTG7XlSCgTFmC9ARjfdSSiWXzQdbWFg+g9IZGe/4WmWeh1996lKq8zP565/t4C8HWiwYoTUuueQSDh8+zLFjxxgaGuKxxx7jpptuitn14pkzuFxEdonIMyIy7rxLRO4WkToRqWttbY3j8JRS8dY7MMzrJzpZXzv+WSwlORk8dvelzC/N4e5f1PHM7lNxHKF1UlJSuP/++7nmmmtYuHAht956K4sXx27JKl7B4A1gpjFmGfB94PfjPdEY85AxZqUxZmVx8ZQP61FKOcjLR9oY8RvWz5/4ez0/K41HPrmaZVV53PvLN/jdGw1xGqG1rr/+eg4dOsRbb73FV77ylZheKy7BwBjTY4zpC378NJAqIkXxuLZSyr42H2wlJz2FiyeRHJ6RkcovPr6Ky+YW8sVf7+SRbdM60EuNIy7BQETKJNhIQ0RWBa/bHo9rK6XsyRjD5oOtrAmWlE5GZloKP/noJVxVW8JXntjDj7cejfEok0e0SksfBdYDRSLSAHwVSAUwxjwA3AJ8WkRGAC9wm7FTRyqlVNwdbA6WlE6QLxhLRqqbBz68gs//6k2+/tR+ur3D/M3VNZMOKJNljLF1s7po/wiNSjAwxtx+nq/fD9wfjWsppRLDaElp2P6CyUpLcXHfbRfiSXPz/T8f4dk9p/mf713EuvPkHiYrIyOD9vZ227axDp1nkJHxzgqs6dJ2FEopS2w+2MKCshzKcqf3Ay3F7eJbtyzj3YtK+ebT+/nIT7dz9YISvnLDQi4ozo5obFVVVTQ0NGDnisbQSWfRosFAKRV3vQPD1B3v5BNrL4jofUSEaxaXsb62mJ+9fJz7/3yEa767hY9eNpvPXV1Drmd6R0KmpqZG7QQxp9DeREqpuBstKZ1ivmA86Slu7rlyLn/50nref1EVP3n5GFd9ezOPbDuBz6/pycnQYKCUirtQSWm0+w0V56Tzb7csY9Nn1zCvOJuvPLGHG763NelaWUyHBgOlVFxNp6R0qpZU5vKrT13KD+64mN6BEe740TY+9Ys6Trb3x+R6iUCDgVIqrqZbUjpVIsINy8p54YtX8qX3zGfr4TY2/PuL/NuzB+gbHInptZ1Ig4FSKq4iKSmdjoxUN599Vw1/+dJ63ru8nB9ufosN33mRzjNDcbm+U2gwUErFVaQlpdNVOiODf7/1Qv7fnRdzumeANxu64np9u9NgoJSKm1BJ6fra+MwKxnLRzDwAmrq8lo3BjjQYKKXi5uUj7VEtKZ2OkpwMUlxCY6cGg3AaDJRScfPioZaYlJROhdsllOVm6MzgHBoMlFJxESopvWJe7EpKJ6siz0OjBoO30WCglIqLQ819nOqOfUnpZFTleWjqGrB6GLaiwUApFRebDwbOLz734HsrVOR5ON0zwIjPb/VQbEODgVIqLjYfbGVBWQ7luR6rh0Jlvgef39DcO2j1UGxDg4FSKuZ6B4bZcbzDFrMCCMwMQMtLw2kwUErF3GhJaZx2HZ9PZTAYaHnpWRoMlFIx9+KhFrLTU1g527qS0nAVeYHdz1pRdJYebqOUw/QMDPNvzxzgZEc/GxaWcs3isri3dpiK0S6lNigpDclMSyE/M1WDQRgNBko5yKtvtfOlX+/kdM8AMwsy+erGvXx1414unpnH9UvLuWZxGdUFmVYP821CJaV/e7U98gUhlfkezRmE0WCglAMMDPv49nMH+fFLx5hTlMVv7rmMi2bmc6Sll2d2n+aZPaf5+lP7+fpT+1lamcu1S8q4bklZxGcBR4OdSkrDVeR6ON5+xuph2IYGA6Vsbk9jN194/E0ONffx3y6dxT9ev4DMtMC37rySHD53dQ6fu7qGE+1neGZPIDB867mDfOu5gywoywkGhnLml2YjInEfv51KSsNV5nt4+UgbxhhL/rvYTVSCgYj8FHgv0GKMWTLG1wW4D7ge6Ac+Zox5IxrXVipRjfj8PLjlKN99/hD5mWk8/NeXTNjtc1ZhFvdcOZd7rpxLU5eXZ/ec5tk9p7nvhcN89/nDXFCUxR2rZ/LxNXPi9sOvb3CEuhMd3LXGfofLV+Z5ODPko8c7Qm5mqtXDsVy0ZgYPA/cDPx/n69cBNcE/q4EfBv9WSo3heNsZvvD4m7xxsosblpXz9ZuXkJ+VNunXV+R5uGvNHO5aM4eW3gH+uLeZ377RwNef2s9VC0qYG6flo5ePtDHss09JabhQeWlDVz+5mbkWj8Z6UUntG2O2AB0TPOVm4Ocm4DUgT0TKo3FtpRKJMYZHtp3g+u9t5UhLH/fddiH3337RlALBuUpyMvjwpbO470MXAbD1UGu0hntemw+22qqkNNzZjWfaowjit8+gEqgP+7wh+Ng7iMjdIlInInWtrfG7aZWyWkvPAHc9vIOvPLGHi2fm89zn13HzhZVRW9KZWZjJ7MJMthxui8r7nY8xhhcPtnDFvELblJSGqxjdeNZv8UjswXb/h4wxDxljVhpjVhYX26v6QKlYeWb3Ka757hZeeaudf75xET+/a1VMEq7r5hfz6lvtDI74ov7e5zrc0kdT94Clp5pNpCg7jbQUF03dOjOA+AWDRqA67POq4GNKJb2vbdrHpx95g6r8TJ76m7V87Io5uFyxSfCurSnGO+zj9ROdMXn/cKGSUju0rB6LiFCp5xqMilcw2Ah8RAIuBbqNMafidG2lbMvnN/z81ePcsLSc333mcuaVxDaxe9ncQlJcwpZDsV8q2nywldpS+5WUhqvM82h/oqCoBAMReRR4FagVkQYR+biI3CMi9wSf8jRwFDgC/Aj4TDSuq5TTner2MuI3rKmJT6uG7PQULp6Vz9bDsc3H9Q2OsON4h21nBSEVeXr8ZUhUSkuNMbef5+sGuDca11IqkdR3BH4QVefHr4XElfOL+dZzB2ntHaQ4Jz0m1xgtKbVpviCkIs9DS+8ggyM+0lPcVg/HUrZLICuVTOo7ApUs1QXxW0pZW1MEBH5gx4qdS0rDhfYanNYksgYDpaxU39mPS86WOcbDkopcCrLS2BKj/QZ2LykNN3qugS4VaTBQykr1Hf2U53ri+kPT5RLWzCtiy+E2/H4T9fff2dBNU/cAVy8sjfp7R1tlvh5yE6LBQCkL1Xd647pEFLK2poi2vkEOnO6N+ntvfLOJNLeLaxaXRf29oy10DoTuQtZgoJSlTnb0M9OC8wfWzQ9U+WyJclWRz2/4w64mrqwtJtdj/+Zv6SluSnLSaezSXcgaDJSyyMCwj9bewbhWEoWUzsigtjQn6iWm24910NI7yE3LK6L6vrFUkefRmQEaDJSyTENnqJLImpPJ1s0vYsexTvqHRqL2nht3NpGZ5maDA/IFIboLOUCDgVIWGd1jYEHOAAKtKYZ8frYdm6jh8OQN+/w8s+cU715UiifNOTX7lfmBYBDYDpW8NBgoZZGTHdbODFbNKSA9xRW1EtOXDrfR1T/Mjcucs0QEUJGbwdCIn/YzQ1YPxVIJFwx8fsNnf/kGv3m9weqhKDWh+o5+MlJdFGfHZhfw+WSkull9QWHUgsHGnU3kelJHk9NOURnM2SR7eWnCBQO3S3jtaDt1x6Mz9VUqVuo7+6nKz7T0/N11NUW81Xom4jVz75CPP+49zXVLykhLcdaPlYq8UHmpBoOEU5WfOToFV8quTnZ4qc63tqNn6Lf4SE8/+8vBFs4M+bjRQVVEIboLOSAhg8HMgkzq9fQiZWPGGBos2mMQrqYkm7IZGWyN8PSzjW82UZyTzqUXFEZpZPGT60klK82twcDqAcRCdUGgbnjE57d6KEqNqds7TO/giGXJ4xARYW1NES8dacM3zdYUPQPD/PlgCzcsLccdo0N5YklEgnsNNBgknOr8THx+wyntRKhsKlRWWmXBhrNzrZ1fTLd3mF0NXdN6/Z/2NjM04nfkElFIqLw0mSVmMAj+tlWveQNlUyctaF09nrXzihBh2qefbdzZRFW+h4tn5kV5ZPGju5ATNBiE1mE1b6Dsqt7i3cfh8rPSWFaZO60+Re19g7x0pI0bl1dYWhUVqco8Dx1nhqK6G9tpEjIYlOdm4HbJ6FRcKbup7+gnLzOVGRn2aOa2tqaYN+u76PYOT+l1T+85jc9vHLfR7FyhiqJknh0kZDBIcbsoz83QmYGyrfpOryUN6sazbn4xPr/h1bemtlS0aWcT80qyWVieE6ORxUeFlpcmZjCAQBJZ9xoou6q3QVlpuItm5pGdnsKWKZSYnur2suN4Bzc5fIkIzh5yk8wVRQkbDGYWZOoykbIlv9/Q2OmlygbJ45BUt4vL5gZaU0y2Ydsfdp7CGBxdRRRSmpOO2yVRDwabdjbxUoR7OOIlKsFARK4VkYMickREvjzG19eLSLeIvBn887+icd2JVBd4aOsbxDvki/WlVBz99c+2850/HrR6GBFp7h1gyOe31TIRBFpTNHR6Od4+uRn1pl1NLK3MZU5RVoxHFnspbhdlMzKi2p/IGMNXN+7ly7/bNe09HPEUcTAQETfwA+A6YBFwu4gsGuOpW40xFwb/fC3S655PqEqjQfMGCcMYw2tHO/jTvmarhxKRs62rbRYMQqefTaI1xbG2M+xq6HbUITbnU5GXEdWcQVP3AB1nhmjo9Eb9RLlYiMbMYBVwxBhz1BgzBDwG3ByF941I6BtN8waJo7N/GO+wj8MtfQwMO3fGF7on7ZQzAJhVmMWswsxJBYNNO5sAeO/y8lgPK26ifcjN7uAmvhSX8MhrJ6P2vrESjWBQCdSHfd4QfOxcl4vILhF5RkQWj/dmInK3iNSJSF1r6/SjaWgKrhvPEkdoPdfnN+w/1WPxaKavvqMfkbPdMu1kbU0Rrx5tZ2hk/FYuxhg27mxi1ewCynPtk/eIVEWeh9PdA1Fb0tnd2E2KS/jY5bP584Fm21cqxSuB/AYw0xizDPg+8PvxnmiMecgYs9IYs7K4ePp90Yuy0/CkuqlP8h7liaQh7P/lnsZuC0cSmfrOfspmZJCeYr/TwNbVFNM/5OP1E53jPufA6V6OtPRx44WJs0QEgWAw4je09EZnr8Huxh5qSnP46OWzMcCvttt7dhCNYNAIVId9XhV8bJQxpscY0xf8+GkgVUSKonDtcYkIVfkeXSZKIKGZQUaqi91ODgYd/bZLHodcNreQFJewdYI17o07m3C7hOuXlMVxZLEXzfJSYwx7GrtZWjmD6oJM1s8v5rEd9QzbuHlmNILBDqBGROaISBpwG7Ax/AkiUibBQmQRWRW8bnsUrj2hQHmpBoNE0dTlJSPVxao5hexqcHIw8NoueRySk5HKxTPzx014GmPYtLOJK+YVUWjRCW2xcvZcg8hnBo1dXjrODLG0MheAO1fPoqV3kBf227f4IeJgYIwZAT4LPAfsBx43xuwVkXtE5J7g024B9ojITuB7wG0mDqdPVxdk0tCpB10nisYuL5V5HpZV5jo2iTw44qO5d8AWDerGs7amiD2NPbT3Db7ja/9V30VDpzehqohCRnchR2FpObSMubQq0LzvqgUlVORm8Mg2+y4VRSVnYIx52hgz3xgz1xjzjeBjDxhjHgh+fL8xZrExZrkx5lJjzCvRuO75VOV76Bscoat/av1WlD01dXmpyPOwpDLXsUnkxk4vxmDbZSI4W2L60pF3bpba+GYTaSku3rO4NN7Dirns9BRyPalRWSYKJY8XlAXadLhdwu2rZrL1cBvH285E/P6xkLA7kEHLSxNNY5eXqnwPS6sCU28nJpFHy0oL7RsMllTmkp+ZyovnlJj6/Iandp/iqtpi2zTYi7aKKJWXhpLHGalniwQ+dEk1bpfwqE0TyQkdDLSVdeIYGPbR1jdERa6HitwMCrLSHJk3CFW32Xlm4HYJV8wrYuvhtrctsW472k5r7yA3LR+rcjwxVEbhxDNjDLsbulhaOeNtj5fMyOA9i0p5vK7elkucCR0Mzh5yo+WlThf6Bq3I8yAiLK3MdWRFUUNHP2kpLkpy7J18XTe/mNbeQQ6c7h19bOPOJrLS3Fy9sMTCkcVWZRR2ITd2eensHx7NF4S7c/UsOvuHeXbP6YiuEQsJHQyy01PIz0zVmUECCPWZD5X/LXVoErm+s5+qPA8um58VvLYmUPkdKjEdGvHzzJ7TvGdx2duWPhJNZb6H3oERegamn2ccTR4HK4nCXT63kNmFmTyy7cS03z9WEjoYQGB2oOWlzheaGYTK/5yaRD7Z0W/bstJw5bke5pdmjx6FufVwK93eYW5MoPYTY6nIi3yvwbnJ43Aul3Dn6lnsON7JwbBZlx1oMFCO0NDlRQTKcgMtHEJJZKctFQX2GNi3rDTcuppith/vwDvkY+POJvIyU1kzb/pdAZygMgrlpbsaut+RPA73gRVVpKW4+KXNZgeJHwzyM2ns8jqihawaX1OXl9KcDFLdgVu2IjeDwqw0djsoidztHabbO2zr5HG4tfOLGRrx8+KhFv60r5nrlpSTlpLYPzIqI5wZhHYeLxtjiSikICuNG5aW87s3Gm115nJi/58lcK7BsM/Q3JO8Z5smgsZO72i+AALtRpY4LIkcmqE6YZkIYPWcAtJSXHzz6QP0D/kSfokIoCg7nTS3i4ZpBoNQ8nhJ1fjBAODO1TPpHRwZ7f5qB4kfDPJ1r0EiaOr2jq7nhjgtiRw6W8NuravHk5HqZvWcAk529FOSk87qOYVWDynmXC6hPC9jtGBhqiZKHodbMSuf2tIc/tNGra0TPhiM7jXQYOBYfr/hVNfAO1o+h5LI+xySRB491MYhy0QQyBsAvHdZBW6bV0BFS0Xu9Pca7GoYP3kcTkS489KZ7G7sZlfw3AOrJXwwCNSlo62sHaytb5Ahn5+qc2YGyxy2E7m+s5+cjBRyM52ze/faJWXMLc7itlXV539ygqjM90w7gby7sZv5EySPw73vokoy09y2Ofgm4YNBWoqL8hkZNOjMwLEawzachSt3WBL5pI1bV4+nuiCTF764nvmlE/+mm0gq8jw09w5Mud302bbVEy8RheRkpHLzhRVs3NlEt9f6/mkJHwwAqgoyNWfgYKFgEJ5ABuclkes7+h2TL0hmVXkejIHT3VPLG0w2eRzujlWz8A77+P1/NZ7/yTGWFMFgZkGm7kJ2sKZxZgbgnCSyMYaGTufsMUhmo62sp5g3CM1QJzszgMB+meVVuTyy7YTlrfaTIhhU52fS3DNo+x8YamyNnV5yMlLG7JTplCRya+8ggyN+x5SVJrNQocJUk8gT7TyeyJ2rZ3GouY+6CY4ajYfkCAYF04v0yh4auwZGNwOdyylJ5JMO22OQzKZ7yM1Uksfh3ru8nJyMFP7zNWt3JCdFMJip5xo4WmPXO/cYhDgliRxapnRaAjkZZaS6KcpOo6l78sFgqsnjcJlpKXzg4iqe2X16zNPl4iUpgkHotzGtKHKmpuBxl2NxShI5tMegKl9zBk5QmeehYQozg+kkj8PduXomQz4/v3m9YVqvj4akCAbF2emkpbh0r4ED9Q2O0O0dHndmAM5IItcHd/EmcvvnRFIxxUNuppM8DldTmsOqOQX8cvtJ/Bb1UUuKYOByCVX5Hk6268zAaZrGKSsNt7TK/knkk1pW6iih4y8nW+Ez3eRxuDtXz+REez8vv/XOs6fjISmCAWh5qVON7jE4pxVFuNBvY3ZOIgfKSjUYOEVlnoeBYT+d/ZPbDDbd5HG4a5eUUZCVZtmO5KgEAxG5VkQOisgREfnyGF8XEfle8Ou7ROTiaFx3Kqrz9VwDJwpVdFTmjf+DNJREtuuZyEMjfk51e6nWfIFjTOWQm0iSx+HSU9x8cGUVf9rfbEmX5YiDgYi4gR8A1wGLgNtFZNE5T7sOqAn+uRv4YaTXnarqAg89AyN0TzLSK3to6vKS4hKKJzgzOJREtuvMoKnLi98EdsIrZwgl+ieTRG7ojCx5HO6OVTPx+Q2Pba+P+L2mKhozg1XAEWPMUWPMEPAYcPM5z7kZ+LkJeA3IE5G4NkcPlfTpUpGzNHZ5KcvNOG/HzGVV9k0i1zusdbWa2swg9EvIRAfaTNaswizWzS/msR0nGZlib6RIRSMYVALhYawh+NhUnwOAiNwtInUiUtfa2hqF4QVUaytrR5qorDRctHci9w2O8PTuU1FpETDaulqDgWPkZ6biSXVPaqNqKHlcG0HyONzn3jWPf7rh3MWV2LNdAtkY85AxZqUxZmVxcfTOWx0NBjozcJSmCXYfh4t2Evn7fz7MZx55gz2NkQeX+s5+Ut1C2Yzxk+DKXkSEiryMSc0MopE8DnfJ7AJuWFZOiju+P56jcbVGILzZeVXwsak+J6ZyPanMyEgZ/S1N2d+Iz8/pnoEJy0pDoplEHhzx8eu6wOafP+1vjvj9Tnb0U5HnSZrDYRJFqLx0IsYYdkcheWwH0QgGO4AaEZkjImnAbcDGc56zEfhIsKroUqDbGHMqCteekmptZe0ozb2D+Pxmwg1nISLC0qroJJGf3XOajjND5HpSeX5f5MGgQfcYOFJV/vk3njV0eunqH2ZpFJLHVos4GBhjRoDPAs8B+4HHjTF7ReQeEbkn+LSngaPAEeBHwGcive506F4DZzlbVjq5ksxo7UR+ZNtJZhVmcs+Vc9l3qmfaRyCG1Hd6qdKeRI5TkeuhrW9owvtpsmceO0FUFqWMMU8bY+YbY+YaY74RfOwBY8wDwY+NMebe4NeXGmPqonHdqaouyKSh02vZdm81NROdYzCWaCSRDzX3sv1YB3esmsm7F5UC8EIES0V9gyN0nBnScwwcKLQ8OdEvA9FOHlvJdgnkWKrO9zA04qfVws6AavLOHnc5ucRr6LezSDqY/nLbSdLcLm5ZUcXc4izmFGXxp/0t036/UPWaLhM5z9ny0vE3gEU7eWylpAoGVdrK2lEau7wUZKWRmZYyqeeX52ZQlJ027Q6m/UMj/PaNBq5bWkZhdjoiwoaFJbz2Vjt9gyPTes9QMNDW1c5TOXri2dg/L0LJ42UJkC+AJAsGM3WvgaM0dXknPSuAyHci/2HnKXoHRvjwpbNGH9uwsJQhn5+th6a35yXUKVf3GDhPWW4GIoHDlcYSSh4vSYB8ASRZMAhFei0vdYbGzsltOAsXSiJ7h6aeRH5k2wnml2azclb+6GMrZuWT60mddolpfUc/WWlu8jPfeWSnsrdUt4vSnPH3GiRS8hiSLBhkpLopnZGuFUUOYIwJzgymFgymm0Te3dDNzoZu7lw9C5Gz+wFS3C7etaCEvxxowTeNwoP6jn6qCzLf9p7KOSrzPeMef7krgZLHkGTBAAJrt5ozsL9u7zBnhnxTnhlM90zkR7adwJPq5n0Xv7NLyoaFpXT2D/PGyakfWF7f2a9LRA5WkecZ9/jLPQmUPIYkDAYzCzL1+EsHOHuOwdSCQdmMqSeRewaGefLNJm5aXsGMjHcu56ybX0SqW6a8Ac0YQ32HV5PHDlaZ5+FU18A7ytETLXkMSRgMqgoyOdUzwNBIfDsCqqkJTc2nukw0nSTy7/+rEe+wjzsvnTnm13MyUrn0gsIp5w3a+obwDvt0j4GDVeZlMOTz03ZOOXqiJY8hCYNBdb4HYybXmtYuBkd8XPN/t/Dkm3Ft52SpyRx3OZ6llbkcau6dVBLZGMMjr51kaWUuy6ryxn3ehoWlHG09w9HWvkmPQ1tXO1/ol5GGc35e7E6w5DEkYTCY6cC9Bq8f7+Rgcy//8cpxq4cSN03dA6SnuCjMSpvya5dW5uI3TCqJ/PqJwH/bD48zKwi5emEJAC9MYQPa6B4DDQaONd4u5ETaeRySdMHAia2sXzwcqHF/42RX0uyRCJWVTqcKZ+kUksiPbDtJTnoKNy6vmPB5VfmZLCjLmdJSUeiUrCo97tKxxjvkZk9jN7VliZM8hiQMBqUzMkh1i6P2Gmw91MYFRVkAbNrVZPFo4qNxGmWlIZNNInecGeKp3ad4/8WVk9rl/O5FpdQd76DzzNCkxnGyvZ/lqdDCAAAajUlEQVSi7PRJ76BW9jMjI5WcjJS3lZcmUtvqcEkXDNwuoTLP45iZQWvvIPtO9fCBFVVcPDOPTTvj3vnbEo2TPOFsLKEk8vl6FP3m9XqGRvzcsXrWhM8L2bCwFL+BzYcmt1QUKCvVWYHTVeZ53rYLORGTx5CEwQACS0VOWW556UhgiWhdTTE3Lq9g/6kejrT0Wjyq2Boc8dHaOzjtmQEEzqM93DJ+EtnvN/xy20kumZ0/6XXfpZW5lOSk8/y+KQQDLSt1vMpzDrlJxOQxaDCwvS2H2ijISmNxxQxuWFaOS2Djm4m9VHQq+FvYVPoSnWvJeZLIr7zVzvH2fu6c5KwAwOUSrl5YwouHWs9bmjzi89PUNaAzgwRQkff2Q252N3aT6hYWlCdO8hiSNRjkZ9LZP0zvwLDVQ5mQ32/YeriNNfOKcLmEkpwMLptbyKZd0Tmo3a4iKSsNOV8S+ZFtJ8jPTOXaJWVTet8NC0vpGxxh27H2CZ93qnsAn99oWWkCqMjz0O0dHu1cG9p5nJ6SOMljSNZgUOCMhnX7T/fQ1jfIuvnFo4/duKyCY21nonJQu11Nd/dxuFASeawzkZt7BvjjvmZuXVk95WqQK+YVkZHqOu9uZG1dnTjCy0uNMexqSLzkMSRpMJjpkPLSLYfaAFhXUzT62HVLykl1Cxt3Ju4GtMYuLyKBFsLTNdFO5Md31OPzG25fNfHegrFkpLpZM6+Y5/e3TDg7C91busfA+SqDy5WNXV4aOr10exMveQxJGgxCv63ZPW+w9XArC8pyKJlx9odibmYqV84v5g+7TiXs8Z1NXV6Ks9MjnoaPlUT2+Q2Pbj/J2poiZgfLdafq3YtKaOzysv/U+In8kx39uF1CeQQBTdlDZV7g50Vjp3c0eZxIPYlCkjIY5GWmkp2eMropyI76h0aoO975tiWikBuXV3Cqe4C6E1PvoukEjV3eiPIFIWMlkf9yoIWm7gHuXD31WUHIuxaUIjLx2cj1HYGDeVLcSfktllCKc9JJcQlNXd7R5HEi7TwOSco7VUSoyvfYuiXFtqMdDPn8rA1bIgrZsLCUjFRXwi4VNXUNRFRWGjJWEvmRbScoyUnn6oWl037f4px0LqzO4/mJgoGWlSYMt0soy82gscvL7obETB5DkgYDCOQN7LxM9OKhVjJSXVwyu+AdX8tKT2HDwlKe3n2aYV9idV/1+01EG87CnZtEru/oZ/OhVm67pJrUCH9j37CwlJ0N3TT3jH0korauTiyVeZ7RZaJETB5DhMFARApE5E8icjj4d/44zzsuIrtF5E0RqYvkmtFSXZBJQ6fXtiWaWw+3snpO4bjVLjcur6DjzBCvvDVxiaPTtJ8ZYmjEH5VgICIsDUsiP7bjJAJ8aBqJ43NtCM4s/nzgnRvQ+odGaOsb1D0GCaQyz8Puxm66vcOjM85EE+nM4MvAC8aYGuCF4OfjucoYc6ExZmWE14yK6nwP3mEfbX2T6zMTT41dXt5qPTPmElHI+tpicjJSEm4DWmiPQTSWiSB0JnIvPQPD/GpHA+9aUBKVQDO/NJvqAs+YJaahXJRWEiWOynwPg8GNhjozGNvNwH8EP/4P4K8ifL+4qbZxK+sthwItKK4cI3kckp7i5trFZfxx72kGhqd++LtdRWOPQbhQEvm+5w/T1jfInZdOfsfxRESEDQtLeelIG/1DI2/7mrauTjyhX04SNXkMkQeDUmNMqHPaaWC8rJwBnheR10Xk7oneUETuFpE6EalrbW2NcHjjC+01aLDhXoOth1spm5HBvJLsCZ934/IKegdH2Hxw8j327a4pysEgNKV/+JXjVOV7WFczfoCdqg0LSxkc8fPS4ba3Pa4bzhJP6H5M1OQxTCIYiMjzIrJnjD83hz/PBBbfx1uAX2OMuRC4DrhXRNaNdz1jzEPGmJXGmJXFxdH7xj1XlU33Goz4Aj9c1s0vOm8v/8vnFlKUnZZQnUwbOr1kp6cwwxOdts+BJHL66CYzt2vq5yOMZ9WcAnIyUt5x4M3JDi+eVDdF2VM/mEfZU2hmkKhLRADn/Y4zxmwY72si0iwi5caYUyJSDoz5K6oxpjH4d4uIPAGsArZMc8xR4UlzU5SdbruWFLsau+kZGBlzf8G5Utwurl9azq921NM3OEJ2uvP75jd1Berzp3OozVgCSeQZbD3cxq0rq6PyniGpbhfra0t44UAzfr/BFQw0odbV0fo3KOtV5XuYX5rNuxdNvyTZ7iJdJtoIfDT48UeBJ899gohkiUhO6GPgPcCeCK8bFdUF9ttrsOVQKyJwxdzxk8fhblxeweCInz/tOx3jkcVHJIfajOfz757P//3QhRTnpEf1fQE2LCyhrW+INxu6Rh+r79A9BokmI9XNHz9/ZUT7U+wu0mDwr8C7ReQwsCH4OSJSISJPB59TCrwkIjuB7cBTxphnI7xuVMwsyLRdf6Kth9tYVpVH/iTP/l0xM5+K3IyEWSpqitIeg3DLqvLOe6zldK2fX4LbJaO7kY0xNHR6NXmsHCeiYGCMaTfGXG2MqTHGbDDGdAQfbzLGXB/8+KgxZnnwz2JjzDeiMfBoqM7P5FT3ACM22bjV7R3mzfqutzWmOx+XS7hxeQVbDrVO+jhGu+ofGqGzfzjqM4NYys1MZdXsgtEDbzr7A62ONRgop0naHcgQWCby+Q2nusfeRRpvrxxpw+c3k8oXhLtxeQUjfsMze5y9VBSqJHLaAfIbFpVysLmXk+39YZVEzvo3KJXkwcBeew22HG4jJz2FC6vzpvS6xRUzuKAoi007nb0BrXH0hDNn/SDdsLAEgOf3N2vrauVYyR0MbFReaoxhy6FWLptbOOW+OSKBpaLXjrWP2yvHCRo7o7vHIF5mFWZRU5LNCweaR3+x0GCgnCapg0F5bgZul9giiXys7QyNXd4pLxGF3Li8AmPgD7ucm0hu6vLidgklMaj6ibUNi0rZdrSDfU09FGSlJUSZr0ouSR0MUtwuKvIybLHXINSCYro7ZOeVZLOofAYbHbxU1NjlpWyGM88A2LCwlBG/4bm9pzVfoBzJed91UTazINMWOYOth9uYXZjJzMLpLy/cdGEFO+u7ONlu/b9nOqLVutoKF1bnUZiVxrDPUKVLRMqBkj4YVOdnWt6faGjEz6tH26e9RBTy3mXlAGza5czZQWj3sRO5XcK7FgQSybrhTDmRBoOCTNr6ht7ReTKe6k500D/kY22ETdSq8jNZMSvfkW2tfX7D6e6BqBx3aZUNwVYFeo6BcqKkDwahmnYr8wZbD7eR4hIum1sY8XvdtLyCg829HDw9/mHtdtTSO8CI3ziurDTc+tpi7r1qLtcuLrN6KEpNWdIHg1ArayvLS7ccauXiWflRqUC5fmk5LsFxew6cWlYaLj3Fzd9fs4DCbOdVQymV9MEgVA9uVXlpW98ge5t6JjzIZiqKc9K5fG4RG3c22fZIz7FE+1AbpdTUJH0wKMxKw5PqtmyZKHQwSjQPXblpeQUnO/rZGTwI3gkao3zcpVJqapI+GIiIpa2stxxupSArjcUVM6L2ntcsKSPVLY5aKmrq8pKXmUqWbtZSyhJJHwwgkDeworzUGMPWw22smVc0ejBKNOR6Urlyfgl/2NWEz++MpaKmrgEqcnVWoJRVNBgQKMms7+iP+xr7/lO9tPYOsnYKLasn66YLK2juGWT7sY6ov3csNHZ6HV1WqpTTaTAgkEQ+M+Sjs384rtfdejjYgiJKyeNwGxaW4El1O2YDWiwOtVFKTZ4GA872no933mDL4VYWlOVQOiP6u24z01JYX1vMn/e32L6qqNs7TO/giAYDpSykwQBG+wHFc6+Bd8jHjmOdMVkiCllfW8zpngEO2HwDWpNWEillOQ0GhJ1rEMck8mvH2hny+WOyRBRy5fxAr5zNB1tjdo1oCG04c2pfIqUSgQYDICs9hYKstLjuNdh6qI30FBeXzC6I2TXKcjNYUJbD5oMtMbtGNDR1BzecaQJZKctoMAiqLsik7ngH7X2DcbnelsOtrL6gkIxUd0yvc9WCEupOdNIzEN/k+FQ0dnlJc7soytI2DkpZJaJgICIfFJG9IuIXkZUTPO9aETkoIkdE5MuRXDNW7rpiNifa+7nmu1t4YX9zTK/V1OXlSEsf62KYLwhZP78Yn9/wcnCnsx01dgZaV0dzr4VSamoinRnsAd4PbBnvCSLiBn4AXAcsAm4XkUURXjfqbr6wko2fu4LinAw+/h91fPm3u+gbjE1b61iWlJ7r4ln55KSn2DpvEDjHQJeIlLJSRMHAGLPfGHPwPE9bBRwxxhw1xgwBjwE3R3LdWFlQNoPf33s5n14/l8fr6rnuvi3sOB79TVtbDrVRNiODmpLsqL/3uVLdLtbUFPHioVbblpg6+YQzpRJFPHIGlUB92OcNwcfGJCJ3i0idiNS1tsb/t9n0FDf//doFPP6pyxCEWx98lX95Zj+DI76ovL/Pb3jpSBtra4oQic+yiJ1LTIdG/LT0DurMQCmLnTcYiMjzIrJnjD8x+e3eGPOQMWalMWZlcXHsl1HGs3J2AU//7Vpuu6SaB188ys33v8yB0z3Tfr/OM0M8vqOej/1sO93e4bgsEYXYucT0dPcAxmjraqWsdt4WkcaYDRFeoxGoDvu8KviY7WWnp/Av71/GhoWl/Pff7uam77/MF98zn0+svQD3JJKdLb0DPLe3mWf3nOK1ox34/IaqfA/3XjWX65bE7zSs8BLTT6+fG7frTsboOQZaVqqUpeLRL3gHUCMicwgEgduAO+Jw3ai5emEpz/1dHv/jid38yzMHeGF/C9+5dfnowTjhGru8PLvnNM/uOUXdiU6MgQuKs7jnygu4bkk5iytmxG15KNz62hJ+vPUovQPD5GSkxv3649Hdx0rZQ0TBQETeB3wfKAaeEpE3jTHXiEgF8GNjzPXGmBER+SzwHOAGfmqM2RvxyOOsMDudBz68gt++0cg/b9zLtd/dwldvXMwHV1ZxsqOfZ/ac5pk9p9lZ3wXAgrIc/vbqGq5fWk5NSbYlASDc+tpiHnjxLV4+0sa1S8otHUu40MygPFd3HytlpYiCgTHmCeCJMR5vAq4P+/xp4OlIrmUHIsItK6q49IICvvTrnfzDb3fx3ecP0dQ9AMCyqlz+4dparltSzpyiLItH+3YrwkpM7RQMmrq8FGWnx3zznVJqYnqs1DRU5Wfyy09cys9eOc6WQ63ctWYO1ywuG3PZyC5S3S6umFfE5oOBElOrZyohgbJSnRUoZTUNBtPkcgkfXzOHj6+ZY/VQJu2qBcU8u/c0B5t7WVAWvWM2I9HY5WVBWY7Vw1Aq6WlvoiRitxJTY0xg97Eed6mU5TQYJJFQielfDtiji2nHmSEGhv1aVqqUDWgwSDLra0t4/UQnvTboYtrUFUi8a1mpUtbTYJBk1tcWM+I3vHzE+i6mjV2Bw4R097FS1tNgkGRW2KiLaWNwZqDBQCnraTBIMueWmFqpsdOLJ9VNXqZ9dkQrlaw0GCShUBfTg83WdjFt6vJSme+xzZ4HpZKZBoMkdGVtoGOq1UtFjXqojVK2ocEgCZXneka7mFqpSQ+1Uco2NBgkqStri6k7bl2J6ZGWXtrPDGkrCqVsQoNBkrqqtiRYYtoe1+t29w/ztU37uPa7W8lJT2F9bUlcr6+UGpv2JkpSoRLTFw+1cG0cDtoZ8fl5dEc9//7Hg3R5h7ntkmq++J5airLTY35tpdT5aTBIUvHsYvrykTa+tmkfB5t7WT2ngP914yIWV+TG7HpKqanTYJDE1tfGtovp8bYzfPPp/fxxXzNV+R5+eOfFXLukTEtJlbIhDQZJLLzENJrBoHdgmPv/fISfvXycFLfw99fU8vE1c/QAG6VsTINBEgsvMb3nyrkRv5/Pb/jN6/V867mDtPUNccuKKv7hmlpKZmjFkFJ2p8EgyV1ZW8xPth6jd2CYnIzpt4XYfqyD/71pL3ubelgxK5+ffuwSllXlRXGkSqlY0tLSJLd+fmQlpsYY7v/zYW598FU6zwzxvdsv4jf3XKaBQCmH0ZlBkls5O5/saZaYGmP4P88d5Ieb3+J9F1XyzfctxZOmeQGlnEiDQZILlJgWTrnE1O83fO0P+3j4lePcuXom/9/NS3C5tEpIKaeKaJlIRD4oIntFxC8iKyd43nER2S0ib4pIXSTXVNG3vraEU90DHGrum9TzfX7Dl3+3i4dfOc4n1szh63+lgUApp4t0ZrAHeD/w4CSee5UxxvrjtdQ7rB8tMW2htixnwucO+/x84fGdbNrZxN9cXcPnN9TovgGlEkBEMwNjzH5jzMFoDUZZ42yJ6cQtrQeGfXz6P99g084mvnzdAr7w7vkaCJRKEPGqJjLA8yLyuojcPdETReRuEakTkbrWVuuPZkwWV9YWU3eig77BkTG/7h3y8cmf1/H8/ma+dvPiqOxLUErZx3mDgYg8LyJ7xvhz8xSus8YYcyFwHXCviKwb74nGmIeMMSuNMSuLi4uncAkVifXzSxj2GV4+8s6VvN6BYT760+28fKSNb92yjI9cNjv+A1RKxdR5cwbGmA2RXsQY0xj8u0VEngBWAVsifV8VPaES080HW7lm8dkS067+IT760+3sberhvtsu4sblFRaOUikVKzFfJhKRLBHJCX0MvIdA4lnZyNkS0xaMMQC09g5y20Ovsf9ULw98eIUGAqUSWKSlpe8TkQbgMuApEXku+HiFiDwdfFop8JKI7AS2A08ZY56N5LoqNsJLTE91e/nQQ69yvP0MP/nYSjYsKrV6eEqpGIqotNQY8wTwxBiPNwHXBz8+CiyP5DoqPkIlpv/52gk2H2qh88wwP79rNavmFFg8MqVUrOkOZDWqPNdDbWkOv3jtBLmeVB75xGqWV2uPIaWSgTaqU2/zgRWVVORm8Njdl2ogUCqJ6MxAvc0n117AJ9deoJvJlEoyGgzU22gQUCo56TKRUkopDQZKKaU0GCillEKDgVJKKTQYKKWUQoOBUkopNBgopZQCJNSh0o5EpBU4AeQC3VN8uZ1fMxM4OcXXTPdadn5NPK+l/6YAvfem/5p4XiuS18wyxkz9MBhjjO3/AA8l2Gta9b+D/ceXoP8mvfcS9//ttP5NoT9OWSbalGCv6ZrGa6Z7LTu/Jp7X0n9TgN57039NPK8Vz38TYPNlokQlInXGmJVWj0MlH7331HicMjNINA9ZPQCVtPTeU2PSmYFSSimdGSillHJ4MBCRKhF5UkQOi8hbInKfiKRN8Py/E5HMeI7RzkTkr0TEiMgCq8cyHhHpO8/XN4tI3NfA9d6bPifcd2DPey+W951jg4EEGu//Dvi9MaYGmA9kA9+Y4GV/B8TtG/J8N5MN3A48Ffx70kTEHZvhOIPeexHT+24aYn7fRVKXauUf4GpgyzmPzQDagSzg28AeYBfwOeBvgCFgN/CXOI2xz+r/ThOMLZvAhr5ZwIHgY+uBLQS+UQ8CDwCu0L8F+A6wE1gTx3H2Bcf1h7DH7gc+Fvx4M7BS7z1n3HtOue/seO/F+r5z8klni4HXwx8wxvSIyEngE8Bs4EJjzIiIFBhjOkTkC8BVxpi2eA1SRLKBJ4F8IBX4J2PMkyIyG3gGeAm4HGgEbjbGeOM0tJuB54wxJ0SkVURWBB9fBSwi8A37LPB+4DcEbrZtxpgvxml8dqb33vTpfTd9Mb3vHLtMdB7rgQeNMSMAxpgOC8cyALzPGHMxcBXwHTl7tmQN8ANjzGICm4E+EMdx3Q48Hvz4cc5O2bcbY44aY3zAo8Ca4OM+4LdxHJ9TrUfvvYnofRcb64nwvnPyzGAfcEv4AyIyg0DvleNWDGgcAnxTRNYBfqASKA1+7Zgx5s3gx68TiOyxH5BIAfAuYKmIGMANGALT9HNrjUOfDwS/Ua0wwtt/ccmwaBwheu9NZzDOu+/AXvdeTO87J88MXgAyReQjMJpc+g7wMPAc8CkRSQl+rSD4ml4gJ87jvBMoBlYYYy4Emjl7Qw2GPc9H/ILzLcAvjDGzjDGzjTHVwDFgLbBKROaIiAv4EIGlBKudABaJSLqI5BFYO7WS3nvT47T7Dux178X0vnNsMDCB7Mn7gA+KyGHgEIFp8f8AfkygM+MuEdkJ3BF82UPAsyLylzgONRdoMcYMi8hVBBJnVrsdeOKcx34bfHwHgSTZfgLfqOc+L26CN/agMaaewJLCnuDf/2XVmEDvvQg44r4De957sb7vdAdyjARvpmaglkADqWygDrgUuC74tD8YY5YEn/8lINsY88/xH22AiKwHvmSMea9VYwgnIsuBHxljVlk9Fidx2r1nt/sOkvPec3LOwO4WA28Fs/iXjfOcJaEPjDHfjsuoHEJE7iFQGvd3Vo/FgfTei0Cy3ns6M4iB8JvJGPNHq8ejkofee2q6NBgopZRybgJZKaVU9GgwiBIR+amItIjInrDHlovIqyKyW0Q2BWuCEZEMEXk0+Ph+EfnHsNesCD5+RES+F7ZJSKkxRfHe+4aI1Nu8r5GKEQ0G0fMwcO05j/0Y+LIxZimBUrm/Dz5+G4EluqXACgL1wbODX/sh8EkCO0RrxnhPpc71MNG59zYRaAuhkpAGgygxxmwBzt0CPp9AAy6AP3F2y/9pICtYAugh0EyqR0TKgRnGmNeCNcU/B/4q5oNXjhaNey/4Pq8ZY07FfsTKjjQYxNZeAo25AD4IVAMYY54FuoFTBDaKfDvYS6QSaAh7fUPwMaWmaqr3nkpyGgxi6y7gMyLyOoEt4UMAIvJhAj3GK4A5wBdF5ALLRqkSkd57akp001kMGWMOAO8BEJH5wA3BL10BPGGMGQZaRORlYCWwFagKe4sqAu2FlZqSadx7Ry0ZqLINnRnEkIiUBP92Af9E4NAOgAMEujciIlkE2gQcCK7X9ojIpcEqoo8Q6Eev1JRM9d6zYozKXjQYRImIPAq8CtSKSIOIfBy4XUQOEfhmawJ+Fnz6g0BasBRwB/AzY8yu4Nc+Q6AS5AjwFoFDSJQaV7TuPRH5PyLSQKAzZoOI/HO8/y3KOroDWSmllM4MlFJKaTBQSimFBgOllFJoMFBKKYUGA6WUUmgwUEophQYDpZRSwP8P4SlnxMf8QcAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "38"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2017-1979"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "38"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2017-1980+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}