fpar_Test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from pyhdf.SD import SD, SDC\n",
    "\n",
    "\n",
    "from netCDF4 import Dataset\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import cv2\n",
    "import os\n",
    "os.getcwd()\n",
    "path1 = 'E:\\\\SIF data\\\\FPAR_2016'\n",
    "path2 = 'E:\\\\SIF data\\\\GOME_2016\\\\GOME_2016\\\\'\n",
    "os.chdir('E:\\\\SIF data\\\\FPAR_2016')\n",
    "# os.chdir('E:\\\\SIF data\\\\OCO_2016')\n",
    "# os.chdir(path2)\n",
    "%matplotlib inline\n",
    "raw_data = SD('FPAR_A2016233.hdf', SDC.READ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from shapely.geometry import Point\n",
    "from shapely.geometry import Polygon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "gnome_test_data = Dataset('08\\\\ret_f_nr5_nsvd12_v26_waves734_nolog.20160801_v27_all.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3678, 10145)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# for i in gnome_test_data.variables:\n",
    "#     print(i)\n",
    "raw_data = SD('FPAR_A2016009.hdf', SDC.READ)\n",
    "raw_data.select('Fpar_500m')[:].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = Polygon([(-88, 40), (-87, 40.5), (-86.5, 41.5), (-87.5, 41)])\n",
    "# p = Polygon([(50, -160), (51, -159), (52, -159), (51, -162)])\n",
    "# ILLINOIS = (36.7281, -91.8524, 42.493703, -87.020001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#find gnome boxes inside illinois\n",
    "\n",
    "lon_corners = gnome_test_data.variables[\"Longitude_corners\"]\n",
    "lat_corners = gnome_test_data.variables[\"Latitude_corners\"]\n",
    "\n",
    "\n",
    "#kd tree?\n",
    "gnome_polygons = []\n",
    "for i in range(len(lon_corners)):\n",
    "    if (lon_corners[i][0] >= ILLINOIS[1] and lon_corners[i][0] <= ILLINOIS[3]) and \\\n",
    "    (lat_corners[i][0] >= ILLINOIS[0] and lat_corners[i][0] <= ILLINOIS[2]):\n",
    "        \n",
    "        #TODO add gnome filter\n",
    "        p1 = (lon_corners[i][0], lat_corners[i][0])\n",
    "        p2 = (lon_corners[i][1], lat_corners[i][1])\n",
    "        p3 = (lon_corners[i][3], lat_corners[i][3])\n",
    "        p4 = (lon_corners[i][2], lat_corners[i][2])\n",
    "               \n",
    "        gnome_box = Polygon([p1, p2, p3, p4])\n",
    "        gnome_polygons.append(gnome_box)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"100.0\" height=\"100.0\" viewBox=\"-88.76392639160156 42.326012725830076 0.7309396362304597 0.5368859863281301\" preserveAspectRatio=\"xMinYMin meet\"><g transform=\"matrix(1,0,0,-1,0,85.18891143798828)\"><path fill-rule=\"evenodd\" fill=\"#66cc99\" stroke=\"#555555\" stroke-width=\"0.014618792724609192\" opacity=\"0.6\" d=\"M -88.22299194335938,42.353084564208984 L -88.06005859375,42.730533599853516 L -88.58002471923828,42.8358268737793 L -88.73685455322266,42.46538162231445 L -88.22299194335938,42.353084564208984 z\" /></g></svg>"
      ],
      "text/plain": [
       "<shapely.geometry.polygon.Polygon at 0x69ca69c9b0>"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gnome_polygons[4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "oco_test_data = Dataset('05\\\\oco2_LtSIF_160504_B7302r_160720173141s.nc4')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([  5.10000000e+01,   9.34000000e+02,   9.27500000e+03,\n",
       "          4.34640000e+04,   5.08100000e+04,   1.57340000e+04,\n",
       "          2.86000000e+03,   2.99000000e+02,   1.10000000e+01,\n",
       "          3.00000000e+00]),\n",
       " array([-2.67721701, -1.97835369, -1.27949038, -0.58062706,  0.11823626,\n",
       "         0.81709957,  1.51596289,  2.2148262 ,  2.91368952,  3.61255283,\n",
       "         4.31141615]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAD8CAYAAACcjGjIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAEUNJREFUeJzt3X+s3XV9x/HnyyJK5hQcV8Paukuy\nZhHJ/NVAF/7ZwEARY9kyEoyTxpE0MZi4xETLTEb8QYJZMjcy50KksTgViT9CI3W1Q4xZIkqZDMHq\nuGNMbiC2WmQYowZ974/76XLWz+k9514K33Ph+UhOzvf7/n6+3/M+zc193e/3+zmnqSokSRr1vKEb\nkCTNHsNBktQxHCRJHcNBktQxHCRJHcNBktQxHCRJHcNBktQxHCRJnZOGbmC1Tj/99Jqfnx+6DUla\nM+6+++4fVdXcNGPXbDjMz89z4MCBoduQpDUjyX9PO9bLSpKkjuEgSeoYDpKkjuEgSeoYDpKkjuEg\nSeoYDpKkjuEgSepMFQ5JHkrynST3JDnQai9Nsj/JA+35tFZPkuuTLCS5N8nrRo6zvY1/IMn2kfrr\n2/EX2r450W9UkjS9lXxC+o+q6kcj6zuB26vquiQ72/p7gYuBTe1xLvAx4NwkLwWuATYDBdydZE9V\nPdbG7ADuBPYCW4EvP6V3pue8+Z23DfK6D113ySCvK51IT+Wy0jZgd1veDVw6Ur+pltwJnJrkDOAi\nYH9VHWmBsB/Y2ra9uKq+UVUF3DRyLEnSAKYNhwK+kuTuJDta7eVV9ShAe35Zq68HHh7Zd7HVlqsv\njqlLkgYy7WWl86rqkSQvA/Yn+d4yY8fdL6hV1PsDLwXTDoBXvOIVy3csSVq1qc4cquqR9nwI+CJw\nDvDDdkmI9nyoDV8ENo7svgF4ZEJ9w5j6uD5uqKrNVbV5bm6qb52VJK3CxHBI8htJfvPoMnAhcB+w\nBzg642g7cGtb3gNc0WYtbQEeb5ed9gEXJjmtzWy6ENjXtj2RZEubpXTFyLEkSQOY5rLSy4Evttml\nJwGfrqp/TnIXcEuSK4EfAJe18XuBNwILwM+AtwNU1ZEkHwTuauM+UFVH2vI7gE8Ap7A0S8mZSpI0\noInhUFUPAq8eU/8xcMGYegFXHedYu4BdY+oHgLOn6FeS9AzwE9KSpI7hIEnqGA6SpI7hIEnqGA6S\npI7hIEnqGA6SpI7hIEnqGA6SpI7hIEnqGA6SpI7hIEnqGA6SpI7hIEnqGA6SpI7hIEnqGA6SpI7h\nIEnqGA6SpI7hIEnqGA6SpI7hIEnqnDR0A3p2m99529AtSFoFzxwkSR3DQZLUMRwkSR3DQZLUMRwk\nSR3DQZLUMRwkSR3DQZLUMRwkSZ2pwyHJuiTfTvKltn5mkm8meSDJZ5Oc3OovaOsLbfv8yDGubvXv\nJ7lopL611RaS7Dxxb0+StBorOXN4F3BwZP3DwEeqahPwGHBlq18JPFZVvwt8pI0jyVnA5cCrgK3A\nP7TAWQd8FLgYOAt4SxsrSRrIVOGQZANwCfDxth7gfOBzbchu4NK2vK2t07Zf0MZvA26uql9U1X8B\nC8A57bFQVQ9W1S+Bm9tYSdJApj1z+FvgPcCv2/pvAT+pqifb+iKwvi2vBx4GaNsfb+P/r37MPser\nS5IGMjEckrwJOFRVd4+WxwytCdtWWh/Xy44kB5IcOHz48DJdS5KeimnOHM4D3pzkIZYu+ZzP0pnE\nqUmOfuX3BuCRtrwIbARo218CHBmtH7PP8eqdqrqhqjZX1ea5ubkpWpckrcbEcKiqq6tqQ1XNs3RD\n+atV9VbgDuBP27DtwK1teU9bp23/alVVq1/eZjOdCWwCvgXcBWxqs59Obq+x54S8O0nSqjyV/+zn\nvcDNST4EfBu4sdVvBD6ZZIGlM4bLAarq/iS3AN8FngSuqqpfASR5J7APWAfsqqr7n0JfkqSnaEXh\nUFVfA77Wlh9kaabRsWN+Dlx2nP2vBa4dU98L7F1JL5Kkp4+fkJYkdQwHSVLHcJAkdQwHSVLHcJAk\ndQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwH\nSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLH\ncJAkdQwHSVJnYjgkeWGSbyX59yT3J3l/q5+Z5JtJHkjy2SQnt/oL2vpC2z4/cqyrW/37SS4aqW9t\ntYUkO0/825QkrcQ0Zw6/AM6vqlcDrwG2JtkCfBj4SFVtAh4DrmzjrwQeq6rfBT7SxpHkLOBy4FXA\nVuAfkqxLsg74KHAxcBbwljZWkjSQieFQS37aVp/fHgWcD3yu1XcDl7blbW2dtv2CJGn1m6vqF1X1\nX8ACcE57LFTVg1X1S+DmNlaSNJCp7jm0v/DvAQ4B+4H/BH5SVU+2IYvA+ra8HngYoG1/HPit0fox\n+xyvLkkayFThUFW/qqrXABtY+kv/leOGteccZ9tK650kO5IcSHLg8OHDkxuXJK3KimYrVdVPgK8B\nW4BTk5zUNm0AHmnLi8BGgLb9JcCR0fox+xyvPu71b6iqzVW1eW5ubiWtS5JWYJrZSnNJTm3LpwBv\nAA4CdwB/2oZtB25ty3vaOm37V6uqWv3yNpvpTGAT8C3gLmBTm/10Mks3rfeciDcnSVqdkyYP4Qxg\nd5tV9Dzglqr6UpLvAjcn+RDwbeDGNv5G4JNJFlg6Y7gcoKruT3IL8F3gSeCqqvoVQJJ3AvuAdcCu\nqrr/hL1DSdKKTQyHqroXeO2Y+oMs3X84tv5z4LLjHOta4Nox9b3A3in6lSQ9A/yEtCSpYzhIkjqG\ngySpYzhIkjqGgySpYzhIkjqGgySpYzhIkjqGgySpYzhIkjqGgySpYzhIkjqGgySpYzhIkjqGgySp\nYzhIkjqGgySpYzhIkjqGgySpYzhIkjqGgySpYzhIkjqGgySpYzhIkjqGgySpYzhIkjqGgySpYzhI\nkjqGgySpYzhIkjqGgySpYzhIkjoTwyHJxiR3JDmY5P4k72r1lybZn+SB9nxaqyfJ9UkWktyb5HUj\nx9rexj+QZPtI/fVJvtP2uT5Jno43K0mazjRnDk8C766qVwJbgKuSnAXsBG6vqk3A7W0d4GJgU3vs\nAD4GS2ECXAOcC5wDXHM0UNqYHSP7bX3qb02StFoTw6GqHq2qf2vLTwAHgfXANmB3G7YbuLQtbwNu\nqiV3AqcmOQO4CNhfVUeq6jFgP7C1bXtxVX2jqgq4aeRYkqQBrOieQ5J54LXAN4GXV9WjsBQgwMva\nsPXAwyO7LbbacvXFMfVxr78jyYEkBw4fPryS1iVJKzB1OCR5EfB54C+q6n+WGzqmVquo98WqG6pq\nc1Vtnpubm9SyJGmVpgqHJM9nKRg+VVVfaOUftktCtOdDrb4IbBzZfQPwyIT6hjF1SdJAppmtFOBG\n4GBV/c3Ipj3A0RlH24FbR+pXtFlLW4DH22WnfcCFSU5rN6IvBPa1bU8k2dJe64qRY0mSBnDSFGPO\nA94GfCfJPa32l8B1wC1JrgR+AFzWtu0F3ggsAD8D3g5QVUeSfBC4q437QFUdacvvAD4BnAJ8uT0k\nSQOZGA5V9a+Mvy8AcMGY8QVcdZxj7QJ2jakfAM6e1Isk6ZnhJ6QlSR3DQZLUMRwkSR3DQZLUMRwk\nSR3DQZLUMRwkSR3DQZLUMRwkSR3DQZLUMRwkSR3DQZLUmeZbWSWtwPzO2wZ53Yeuu2SQ19Wzk2cO\nkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO\n4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6kwMhyS7khxKct9I7aVJ9id5oD2f1upJcn2ShST3JnndyD7b\n2/gHkmwfqb8+yXfaPtcnyYl+k5KklZnmzOETwNZjajuB26tqE3B7Wwe4GNjUHjuAj8FSmADXAOcC\n5wDXHA2UNmbHyH7HvpYk6Rk2MRyq6uvAkWPK24DdbXk3cOlI/aZacidwapIzgIuA/VV1pKoeA/YD\nW9u2F1fVN6qqgJtGjiVJGshq7zm8vKoeBWjPL2v19cDDI+MWW225+uKYuiRpQCf6hvS4+wW1ivr4\ngyc7khxIcuDw4cOrbFGSNMlqw+GH7ZIQ7flQqy8CG0fGbQAemVDfMKY+VlXdUFWbq2rz3NzcKluX\nJE2y2nDYAxydcbQduHWkfkWbtbQFeLxddtoHXJjktHYj+kJgX9v2RJItbZbSFSPHkiQN5KRJA5J8\nBvhD4PQkiyzNOroOuCXJlcAPgMva8L3AG4EF4GfA2wGq6kiSDwJ3tXEfqKqjN7nfwdKMqFOAL7eH\nTrD5nbcN3YKkNWRiOFTVW46z6YIxYwu46jjH2QXsGlM/AJw9qQ9J0jPHT0hLkjqGgySpYzhIkjqG\ngySpYzhIkjqGgySpYzhIkjqGgySpYzhIkjqGgySpYzhIkjqGgySpYzhIkjqGgySpYzhIkjqGgySp\nYzhIkjqGgySpYzhIkjqGgySpYzhIkjqGgySpYzhIkjqGgySpc9LQDUg6MeZ33jbYaz903SWDvbae\nHp45SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6M/M5hyRbgb8D1gEfr6rrBm7phBtyHrok\nrcRMnDkkWQd8FLgYOAt4S5Kzhu1Kkp67ZiIcgHOAhap6sKp+CdwMbBu4J0l6zpqVy0rrgYdH1heB\ncwfqRdIKDXXJ1K/tePrMSjhkTK26QckOYEdb/WmSHwM/ejobO4FOZ+30Cmur37XUK6ytfme613y4\nK810v2M80/3+zrQDZyUcFoGNI+sbgEeOHVRVNwA3HF1PcqCqNj/97T11a6lXWFv9rqVeYW31u5Z6\nBfs9kWblnsNdwKYkZyY5Gbgc2DNwT5L0nDUTZw5V9WSSdwL7WJrKuquq7h+4LUl6zpqJcACoqr3A\n3hXudsPkITNjLfUKa6vftdQrrK1+11KvYL8nTKq6+76SpOe4WbnnIEmaIWs+HJJ8MMm9Se5J8pUk\nvz10T8eT5K+TfK/1+8Ukpw7d03KSXJbk/iS/TjKTMyqSbE3y/SQLSXYO3c9ykuxKcijJfUP3MkmS\njUnuSHKw/Qy8a+ielpPkhUm+leTfW7/vH7qnSZKsS/LtJF8aupdx1nw4AH9dVb9fVa8BvgT81dAN\nLWM/cHZV/T7wH8DVA/czyX3AnwBfH7qRcdbg1658Atg6dBNTehJ4d1W9EtgCXDXj/7a/AM6vqlcD\nrwG2JtkycE+TvAs4OHQTx7Pmw6Gq/mdk9TcY8+G5WVFVX6mqJ9vqnSx9nmNmVdXBqvr+0H0sY019\n7UpVfR04MnQf06iqR6vq39ryEyz9Els/bFfHV0t+2laf3x4z+7sgyQbgEuDjQ/dyPGs+HACSXJvk\nYeCtzPaZw6g/B748dBNr3LivXZnZX2BrVZJ54LXAN4ftZHntMs09wCFgf1XNcr9/C7wH+PXQjRzP\nmgiHJP+S5L4xj20AVfW+qtoIfAp45yz32sa8j6XT9k8N1+n/9TKx3xk21deuaPWSvAj4PPAXx5yl\nz5yq+lW7vLwBOCfJ2UP3NE6SNwGHquruoXtZzsx8zmE5VfWGKYd+GrgNuOZpbGdZk3pNsh14E3BB\nzcA84hX8286iqb52RauT5PksBcOnquoLQ/czrar6SZKvsXR/ZxZv/p8HvDnJG4EXAi9O8k9V9WcD\n9/X/rIkzh+Uk2TSy+mbge0P1Mkn7D43eC7y5qn42dD/PAn7tytMkSYAbgYNV9TdD9zNJkrmjs/+S\nnAK8gRn9XVBVV1fVhqqaZ+ln9quzFgzwLAgH4Lp2GeRe4EKWZgDMqr8HfhPY36be/uPQDS0nyR8n\nWQT+ALgtyb6hexrVbu4f/dqVg8Ats/y1K0k+A3wD+L0ki0muHLqnZZwHvA04v/2s3tP+0p1VZwB3\ntN8Dd7F0z2Emp4iuFX5CWpLUeTacOUiSTjDDQZLUMRwkSR3DQZLUMRwkSR3DQZLUMRwkSR3DQZLU\n+V89hfE9AdBN3gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x69860dfb00>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(oco_sifs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#find oco boxes inside illinois\n",
    "\n",
    "\n",
    "area = []\n",
    "desired_sif = []\n",
    "\n",
    "oco_lat_corners = oco_test_data.variables['footprint_vertex_latitude'][:]\n",
    "oco_lon_corners = oco_test_data.variables['footprint_vertex_longitude'][:]\n",
    "oco_sifs = oco_test_data.variables['SIF_757nm'][:]\n",
    "\n",
    "\n",
    "oco_illinois = []\n",
    "oco_illinois_sifs = []\n",
    "for i in range(len(oco_lon_corners)):\n",
    "    if (oco_lon_corners[i][0] >= ILLINOIS[1] and oco_lon_corners[i][0] <= ILLINOIS[3]) and \\\n",
    "    (oco_lat_corners[i][0] >= ILLINOIS[0] and oco_lat_corners[i][0] <= ILLINOIS[2]):\n",
    "        \n",
    "        \n",
    "        #add filter\n",
    "        #add cdl data filter\n",
    "        oco_illinois.append([oco_lon_corners[i], oco_lat_corners[i]])\n",
    "        oco_illinois_sifs.append(oco_sifs[i])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "area = []\n",
    "desired_sif = []\n",
    "chosen_area = None\n",
    "\n",
    "for j in range(len(gnome_polygons)):\n",
    "    g_area = gnome_polygons[j]\n",
    "    desired_sif = []\n",
    "    area = []\n",
    "    \n",
    "    for i in range(len(oco_illinois)):\n",
    "        p = Point(oco_illinois[i][0][0], oco_illinois[i][1][0])\n",
    "        if g_area.contains(p):\n",
    "            desired_sif.append(oco_illinois_sifs[i])\n",
    "            area.append([oco_illinois[i][0], oco_illinois[i][1]])\n",
    "\n",
    "    if area != []:\n",
    "        chosen_area = g_area\n",
    "        break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "verifying\n",
      "202 202\n",
      "85 85\n",
      "POLYGON ((-87.12796783447266 42.12539672851563, -86.95050048828125 42.52025985717773, -88.55459594726563 42.90095138549805, -88.71170043945313 42.53056335449219, -87.12796783447266 42.12539672851563))\n"
     ]
    },
    {
     "data": {
      "image/svg+xml": [
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"100.0\" height=\"100.0\" viewBox=\"-88.7821484375 42.05494873046875 1.9020959472656216 0.9164506530761685\" preserveAspectRatio=\"xMinYMin meet\"><g transform=\"matrix(1,0,0,-1,0,85.02634811401367)\"><path fill-rule=\"evenodd\" fill=\"#66cc99\" stroke=\"#555555\" stroke-width=\"0.03804191894531243\" opacity=\"0.6\" d=\"M -87.12796783447266,42.125396728515625 L -86.95050048828125,42.520259857177734 L -88.55459594726562,42.90095138549805 L -88.71170043945312,42.53056335449219 L -87.12796783447266,42.125396728515625 z\" /></g></svg>"
      ],
      "text/plain": [
       "<shapely.geometry.polygon.Polygon at 0x69ca69c8d0>"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAEv1JREFUeJzt3X2sXHd95/H3p3loVUiXpL4JIYlx\nH0zUFC2BvTLQCAgPSYNBhFa0ddQHs03rlgepaLur9ZaKrOCftBWg7QY1eImVUEEa0RKwGofETUGB\nKobcBOeJAA7ZtDG2YtNAQkS7rdtv/5jj7eQy13c8Z3zHzu/9kkZzzu/8zvy+nnvzuSe/OXNOqgpJ\nUjt+YNYFSJJWlsEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JasyJsy5glFWrVtWa\nNWtmXYYkHTfuuuuub1XV3Dh9j8ngX7NmDQsLC7MuQ5KOG0n+dty+TvVIUmMMfklqjMEvSY0x+CWp\nMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1Jjjslv7krHqjWbb5rJuI9c+YaZjKtnJo/4JakxBr8kNcbg\nl6TGGPyS1BiDX5Ias2zwJzknyWeTPJjkgSS/07WflmRHkt3d86lL7L+x67M7ycZp/wMkSUdmnCP+\ng8DvVtVPAS8D3pHkPGAzcFtVrQVu69afJslpwBXAS4F1wBVL/YGQJK2MZYO/qvZV1d3d8neBB4Gz\ngEuB67pu1wFvHrH7zwI7qurxqvo2sAO4ZBqFS5Imc0Rz/EnWAC8GvgicUVX7YPDHATh9xC5nAY8O\nre/p2iRJMzJ28Cd5NvAXwLuq6slxdxvRVku8/qYkC0kWDhw4MG5ZkqQjNFbwJzmJQeh/rKo+2TU/\nluTMbvuZwP4Ru+4BzhlaPxvYO2qMqtpSVfNVNT83N9aN4iVJExjnrJ4A1wAPVtUHhjZtAw6dpbMR\n+PSI3W8BLk5yaveh7sVdmyRpRsY54r8A+FXgNUl2dY/1wJXARUl2Axd16ySZT/IRgKp6HHgfcGf3\neG/XJkmakWWvzllVX2D0XD3Aa0f0XwB+Y2h9K7B10gIlSdPlN3clqTEGvyQ1xuCXpMYY/JLUGINf\nkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY1Z9kYs\nSbYCbwT2V9ULu7YbgHO7Ls8BvlNV54/Y9xHgu8C/AAeran5KdUuSJrRs8APXAlcBHz3UUFW/dGg5\nyfuBJw6z/6ur6luTFihJmq5xbr14e5I1o7Z1N2L/ReA10y1LknS09J3jfwXwWFXtXmJ7AbcmuSvJ\npsO9UJJNSRaSLBw4cKBnWZKkpfQN/suA6w+z/YKqegnweuAdSV65VMeq2lJV81U1Pzc317MsSdJS\nJg7+JCcCPw/csFSfqtrbPe8HbgTWTTqeJGk6+hzxvw74alXtGbUxybOSnHJoGbgYuL/HeJKkKVg2\n+JNcD9wBnJtkT5LLu00bWDTNk+R5SbZ3q2cAX0hyD/Al4Kaq+sz0SpckTWKcs3ouW6L9rSPa9gLr\nu+WHgRf1rE8SsGbzTTMb+5Er3zCzsXV0+M1dSWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5Ia\nY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGjPOjVi2Jtmf5P6htv+Z5JtJdnWP\n9Uvse0mSryV5KMnmaRYuSZrMOEf81wKXjGj/YFWd3z22L96Y5ATgQwxutH4ecFmS8/oUK0nqb9ng\nr6rbgccneO11wENV9XBV/RPwZ8ClE7yOJGmK+szxvzPJvd1U0Kkjtp8FPDq0vqdrkyTN0KTB/yfA\nTwDnA/uA94/okxFttdQLJtmUZCHJwoEDByYsS5K0nImCv6oeq6p/qap/Bf4Pg2mdxfYA5wytnw3s\nPcxrbqmq+aqan5ubm6QsSdIYJgr+JGcOrf4ccP+IbncCa5P8WJKTgQ3AtknGkyRNz4nLdUhyPXAh\nsCrJHuAK4MIk5zOYunkE+K2u7/OAj1TV+qo6mOSdwC3ACcDWqnrgqPwrJEljWzb4q+qyEc3XLNF3\nL7B+aH078H2nekqSZsdv7kpSYwx+SWqMwS9JjTH4JakxBr8kNWbZs3oktW3N5ptmMu4jV75hJuO2\nwCN+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY1ZNviTbE2yP8n9\nQ21/lOSrSe5NcmOS5yyx7yNJ7kuyK8nCNAuXJE1mnCP+a4FLFrXtAF5YVf8R+DrwPw6z/6ur6vyq\nmp+sREnSNC0b/FV1O/D4orZbq+pgt7oTOPso1CZJOgqmMcf/68DNS2wr4NYkdyXZNIWxJEk99bos\nc5J3AweBjy3R5YKq2pvkdGBHkq92/wcx6rU2AZsAVq9e3acsSdJhTHzEn2Qj8Ebgl6uqRvWpqr3d\n837gRmDdUq9XVVuqar6q5ufm5iYtS5K0jImCP8klwH8H3lRV31uiz7OSnHJoGbgYuH9UX0nSyhnn\ndM7rgTuAc5PsSXI5cBVwCoPpm11Jru76Pi/J9m7XM4AvJLkH+BJwU1V95qj8KyRJY1t2jr+qLhvR\nfM0SffcC67vlh4EX9apOkjR1fnNXkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiD\nX5Ia0+vqnNKsrNl806xLkI5bHvFLUmMMfklqjMEvSY0x+CWpMQa/JDXG4JekxowV/Em2Jtmf5P6h\nttOS7Eiyu3s+dYl9N3Z9dnf36ZUkzdC4R/zXApcsatsM3FZVa4HbuvWnSXIacAXwUgY3Wr9iqT8Q\nkqSVMVbwV9XtwOOLmi8FruuWrwPePGLXnwV2VNXjVfVtYAff/wdEkrSC+szxn1FV+wC659NH9DkL\neHRofU/X9n2SbEqykGThwIEDPcqSJB3O0f5wNyPaalTHqtpSVfNVNT83N3eUy5KkdvUJ/seSnAnQ\nPe8f0WcPcM7Q+tnA3h5jSpJ66hP824BDZ+lsBD49os8twMVJTu0+1L24a5Mkzci4p3NeD9wBnJtk\nT5LLgSuBi5LsBi7q1kkyn+QjAFX1OPA+4M7u8d6uTZI0I2NdlrmqLlti02tH9F0AfmNofSuwdaLq\nJElT5zd3JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8\nktQYg1+SGmPwS1JjDH5JaszEwZ/k3CS7hh5PJnnXoj4XJnliqM97+pcsSepjrBuxjFJVXwPOB0hy\nAvBN4MYRXT9fVW+cdBxJ0nRNa6rntcA3qupvp/R6kqSjZFrBvwG4foltL09yT5Kbk/z0lMaTJE2o\nd/AnORl4E/CJEZvvBp5fVS8C/jfwqcO8zqYkC0kWDhw40LcsSdISpnHE/3rg7qp6bPGGqnqyqp7q\nlrcDJyVZNepFqmpLVc1X1fzc3NwUypIkjTKN4L+MJaZ5kjw3Sbrldd14fz+FMSVJE5r4rB6AJD8M\nXAT81lDbbwNU1dXAW4C3JTkI/AOwoaqqz5iSpH56BX9VfQ/40UVtVw8tXwVc1WcMSdJ09Qp+tW3N\n5ptmXYKkCXjJBklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfgl\nqTEGvyQ1xuCXpMYY/JLUGINfkhozjZutP5LkviS7kiyM2J4kf5zkoST3JnlJ3zElSZOb1o1YXl1V\n31pi2+uBtd3jpcCfdM+SpBlYiameS4GP1sBO4DlJzlyBcSVJI0zjiL+AW5MU8OGq2rJo+1nAo0Pr\ne7q2fcOdkmwCNgGsXr16CmW1w1sg6plolr/Xj1z5hpmNvRKmccR/QVW9hMGUzjuSvHLR9ozYp76v\noWpLVc1X1fzc3NwUypIkjdI7+Ktqb/e8H7gRWLeoyx7gnKH1s4G9fceVJE2mV/AneVaSUw4tAxcD\n9y/qtg34te7snpcBT1TVPiRJM9F3jv8M4MYkh17r41X1mSS/DVBVVwPbgfXAQ8D3gP/cc0xJUg+9\ngr+qHgZeNKL96qHlAt7RZxxJ0vT4zV1JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtS\nYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1Jhp3HNXkp5RZnW/35W61+/ER/xJzkny2SQPJnkg\nye+M6HNhkieS7Ooe7+lXriSprz5H/AeB362qu7vbL96VZEdVfWVRv89X1Rt7jCNJmqKJj/iral9V\n3d0tfxd4EDhrWoVJko6OqXy4m2QN8GLgiyM2vzzJPUluTvLT0xhPkjS53h/uJnk28BfAu6rqyUWb\n7waeX1VPJVkPfApYu8TrbAI2AaxevbpvWZKkJfQ64k9yEoPQ/1hVfXLx9qp6sqqe6pa3AyclWTXq\ntapqS1XNV9X83Nxcn7IkSYfR56yeANcAD1bVB5bo89yuH0nWdeP9/aRjSpL66zPVcwHwq8B9SXZ1\nbb8HrAaoqquBtwBvS3IQ+AdgQ1VVjzElST1NHPxV9QUgy/S5Crhq0jEkSdPnJRskqTEGvyQ1xuCX\npMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNeYZd8/dWd0rU5KOFx7xS1JjDH5JaozB\nL0mNMfglqTEGvyQ1pu89dy9J8rUkDyXZPGL7Dya5odv+xSRr+ownSeqvzz13TwA+BLweOA+4LMl5\ni7pdDny7qn4S+CDwB5OOJ0majj5H/OuAh6rq4ar6J+DPgEsX9bkUuK5b/nPgtYduvi5Jmo0+wX8W\n8OjQ+p6ubWSfqjoIPAH8aI8xJUk99fnm7qgj95qgz6BjsgnY1K0+leRrPWrrYxXwrRmNfaSOl1qt\nc7qsc7qOmTpz+Mnw5ep8/rjj9An+PcA5Q+tnA3uX6LMnyYnAfwAeH/ViVbUF2NKjnqlIslBV87Ou\nYxzHS63WOV3WOV0t1tlnqudOYG2SH0tyMrAB2LaozzZgY7f8FuCvq2rkEb8kaWVMfMRfVQeTvBO4\nBTgB2FpVDyR5L7BQVduAa4A/TfIQgyP9DdMoWpI0uV5X56yq7cD2RW3vGVr+R+AX+owxAzOfbjoC\nx0ut1jld1jldzdUZZ14kqS1eskGSGtN88Cc5LcmOJLu751OX6PeHSR5I8mCSP57FF9GOoNbVSW7t\nav3KSl8qY9w6u74/kuSbSa5ayRq7sZetM8n5Se7ofvb3JvmlFazvuLgkyhh1/pfu9/DeJLclGfu0\nw5Wsc6jfW5JUkpmc6TNOnUl+sXtPH0jy8SMepKqafgB/CGzuljcDfzCiz88Af8PgQ+wTgDuAC4/F\nWrttnwMu6pafDfzwsVhnt/1/AR8HrjoW30/gBcDabvl5wD7gOStQ2wnAN4AfB04G7gHOW9Tn7cDV\n3fIG4IYZvIfj1PnqQ7+DwNuO1Tq7fqcAtwM7gfljsU5gLfBl4NRu/fQjHaf5I36eflmJ64A3j+hT\nwA8x+EH8IHAS8NiKVPd0y9baXS/pxKraAVBVT1XV91auRGC895Qk/wk4A7h1hepabNk6q+rrVbW7\nW94L7AfmVqC24+WSKMvWWVWfHfod3MngOz8rbZz3E+B9DA4I/nElixsyTp2/CXyoqr4NUFX7j3QQ\ngx/OqKp9AN3z6Ys7VNUdwGcZHO3tA26pqgdXtMqBZWtlcIT6nSSfTPLlJH/UXVBvJS1bZ5IfAN4P\n/LcVrm3YOO/n/5dkHYM//t9YgdqOl0uijFPnsMuBm49qRaMtW2eSFwPnVNVfrmRhi4zzfr4AeEGS\nv0myM8klRzrIM+5m66Mk+SvguSM2vXvM/X8S+Cn+/UhlR5JXVtXtUypxeKxetTL4mb4CeDHwd8AN\nwFsZfKdiaqZQ59uB7VX16NE8SJ1CnYde50zgT4GNVfWv06htuSFHtE18SZSj6Eguy/IrwDzwqqNa\n0WiHrbM7EPkgg/9WZmmc9/NEBtM9FzLIpM8neWFVfWfcQZoI/qp63VLbkjyW5Myq2tf9xz3qf5t+\nDthZVU91+9wMvIzBXOCxVuse4MtV9XC3z6e6Wqca/FOo8+XAK5K8ncHnECcneaqqlvzQbUZ1kuRH\ngJuA36+qndOs7zCmekmUo2icOknyOgZ/bF9VVf9vhWobtlydpwAvBD7XHYg8F9iW5E1VtbBiVY7/\nc99ZVf8M/N8Mrmu2lsHVFMbiVM/TLyuxEfj0iD5/B7wqyYlJTmJwxDKLqZ5xar0TODXJoXno1wBf\nWYHahi1bZ1X9clWtrqo1wH8FPjrt0B/DsnV2lyO5kUF9n1jB2o6XS6IsW2c3hfJh4E2TzEdPyWHr\nrKonqmpVVa3pfid3Mqh3JUN/2To7n2LwgTlJVjGY+nn4iEZZ6U+tj7UHgznR24Dd3fNpXfs88JH6\n90/aP8wg7L8CfOBYrbVbvwi4F7gPuBY4+Visc6j/W5nNWT3j/Ox/BfhnYNfQ4/wVqm898HUGnym8\nu2t7L4NAgsEJB58AHgK+BPz4jH4vl6vzrxicDHHo/dt2LNa5qO/nmMFZPWO+nwE+0GXRfcCGIx3D\nb+5KUmOc6pGkxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ15t8AYQjJNH4XH/gAAAAA\nSUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x698a84ec88>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "print('verifying')\n",
    "\n",
    "print(len(oco_illinois), len(oco_illinois_sifs))\n",
    "\n",
    "print(len(desired_sif), len(area))\n",
    "\n",
    "plt.hist(desired_sif)\n",
    "\n",
    "print(chosen_area)\n",
    "chosen_area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[array([-87.13839722, -87.14543152, -87.1386261 , -87.13159943], dtype=float32), array([ 42.13482285,  42.15303802,  42.14352417,  42.12530899], dtype=float32)]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAD8CAYAAABpcuN4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xd8FHX+x/HXJ400CCWhBlIA6TWh\nJYDYDk+UIiigVClKsZ7lp2c5r1jPsyBIla4gIOUUsWAhJLREehVSCD0QEkogpHx/f2TxIiQhCUlm\nk/08H499uDsz35k3cXc/OzPf+Y4YY1BKKaXy42R1AKWUUvZNC4VSSqkCaaFQSilVIC0USimlCqSF\nQimlVIG0UCillCqQFgqllFIF0kKhlFKqQFoolFJKFcjF6gAlwdfX1wQGBlodQymlypWYmJjTxhi/\nGy1XIQpFYGAg0dHRVsdQSqlyRUQSCrOcHnpSSilVIC0USimlCqSFQimlVIG0UCillCqQFgqllFIF\n0kKhlFKqQFoolFJKFcihC8WlK1nM35jAsZRLVkdRSim75dCFIiH5Iq+u3EX3d37iqUVb2XU01epI\nSilldxy6UDStXYUZQ0Op5OLEim3HuHfSeh6euZGf95/CGGN1PKWUsgsOXSgA7mxei6XjwqhX1QOA\nyINnGDF7C3d/EMGS6ETSM7MsTqiUUtZy+EIB0KxOFVZMCKddg6oANPTzQgSeW7qDbm//xJSfD5Ka\nlmFxSqWUsoYWChu/ypX4fExn+rSty6GkizSrU4WZw0JpUrsy76zZT5e31vL6f3eTmJxmdVSllCpT\nFWL02JLi7urMBwPb0sjPm/e+P8Dh5DSmDQ3h1Ll0ZkbEMn9DAnOj4rmnVR3Gdg+mtX9VqyMrpVSp\nk4pw0jY0NNSU9DDjX+84zl+WbMPXuxKzhnegSe3KHE+9xJyoeD7beJjz6Zl0CqrOmG7B3N60Jk5O\nUqLbV0qp0iYiMcaY0Bsup4UifzuOpDBmXjQX07OYNLgdtzWtCcD5yxks3pLIp+vjOJZ6mWA/L8Z0\nC6Zfu3q4uzqXeA6llCoNWihKyInUy4yet4U9x87x0j3NGNU1CJGcvYeMrGxW7zzO9HWx7D52Dl9v\nN4Z1CWRI5wCqe7mVSh6llCopWihKUNqVTJ5ZvJ01u08wuGN9Xu/dEjeX//UDMMawIfYMM9bF8tP+\nJNxdnXggpD6jugYR6OtVarmUUupmaKEoYdnZhve+38/knw7RObg6U4eEUNXz+r2GAyfPMzMilhVb\nj5GRnU3P5rUZ0z2YkIBqpZpPKaWKSgtFKVm+9QgvLN1J3aruzBrRgYZ+3nkud+rcZeZuiGfBxsOk\nXsogJKAaY7oFc1fzWjjriW+llB3QQlGKYhKSGTsvhoysbKY8HELXxr75LnsxPZMl0YnMiowjMfkS\ngTU8GdUtmAHt/fFw0xPfSinraKEoZYnJaYyeG83BpAu83rsFQzoHFLh8ZlY23+4+yfR1h9h+JJVq\nnq4M7RLIsC4B+HpXKqPUSin1P1ooysD5yxk8uWgbP+47xYiwQF7u1QwX54IvdjfGsCX+LDMiYvlh\n70lcnZ3o396f0d2C8j2MpZRSpUELRRnJyja8uXovM9fHcestfkx6qB1V3F0L1fZQ0gVmrY9jacwR\nrmRmc2ezWozpFkTHoOq/d8FVSqnSooWijC3afJiXV+wi0NeLT4d3oEENz0K3PX0hnfkbEpi3IZ6z\naRm08fdhTPdg7m5R+4Z7KEopVVxaKCwQdeg04xb8ipPAtKGhdAyqXqT2l65ksezXI8yMiCX+TBr+\n1TwY1TWIB0Pr41VJh+VSSpWswhaKQv9cFRFnEdkqIl/ZXi8Ukf0isktEPhWRPI+3iMgaEUm52i7X\n9DkiEici22yPtrbpIiIfichBEdkhIu0Lm9FqYQ19WTEhnGqebjw8cyNLohOL1N7DzZkhnQNY+5ce\nTBsaQu0q7rz+3z10eXMt76zZx6lzl0spuVJK5a8oxzWeBPbmer0QaAq0AjyA0fm0excYms+854wx\nbW2PbbZpfwYa2x5jgU+KkNFyQb5eLB8fTqegGjy3dAdvfrOX7Oyi7bU5Owk9W9Rm6bgwlo0LI7yR\nL5/8coiub//Ec0u2c+Dk+VJKr5RS1ytUoRARf6AXMPPqNGPMamMDbAb882prjFkLFOWbrQ8wz7bq\njUBVEalThPaW8/F0ZfbIDgzp3IBpv8Ty6IIYLqZnFmtdIQHV+GRICD8/24NBHevz3x3H+NP76xgx\nezNRB0/rLVuVUqWusHsUHwDPA9nXzrAdchoKrCnG9v9lO7z0vohcvZigHpD7mM0R27RyxdXZiX/0\nacnrvVuwdu9JBkzdwLGUS8VeX0ANL/7epyUb/u8O/nLXLew6mspDMzdx76T1rNx2lIys6/7XKKVU\nibhhoRCRe4FTxpiYfBaZAqwzxkQUcdsvknPoqgNQHXjh6ibzWPa6n80iMlZEokUkOikpqYibLhsi\nwvCwQGaP7MiR5DR6fxzJ1sNnb2qd1bzcePyOxqx/4Xbeur8VlzOyeHLRNm595ydmRsRy/rLeslUp\nVbIKs0cRDvQWkXhgEXC7iCwAEJHXAD/gmaJu2Bhz3HZ4KR2YDXS0zToC1M+1qD9wLI/2040xocaY\nUD8/v6JuvkzdeosfX44Pw9PNmYHTN7Jq+3X/nCJzd3VmUMcGfP/0rcwaHkr96p788+u9hL35I2+u\n3svx1OLvvSilVG43LBTGmBeNMf7GmEBgEPCjMWaIiIwGegKDjTFFPu5x9byD5FxZ1hfYZZu1Chhm\n6/3UGUg1xhwv6vrtTeNalVkxIZy2/lV54vOt/Of7AyVyfsHJSbijWS0WP9qFVRPD6dG0JjPXx9Ht\n7Z94evE2dh9LLYH0SilHVqTrKESkB/CsMeZeEckEEvjfieovjTF/F5FQ4DFjzGhbmwhyDjF5A2eA\nUcaYb0XkR3L2RgTYZmtzwVY4PgbuBtKAkcaYAi+SsJfrKAojPTOLvy7fxdKYI/RqXYf3HmhT4nfF\nS0xOY3ZkPIu2HCbtShZdG/kypnsw3Rv76hXfSqnf6QV3dswYw/R1sby1Zh+t6/kwY1goNau4l/h2\nUtMy+GzzYWZHxnHqfDpNa1dmdLdgerep+4cbLymlHJMWinLgu90neGrxNnw8XJkxLJSW9XxKZTtX\nMrNZtf0YM9bFsv/keWpVqcSIsCAe6tQAH4/CjUullKp4tFCUE7uPpTJmbjRn0zL4YFBberaoXWrb\nMsaw7rfTzFgXy/qDp/Fyc2ZghwaMDA+kfvXCj02llKoYtFCUI6fOX2bsvBi2H0nh+Z5NeezW4FI/\nl7D7WCozI+L47/ZjGOCeVnUY0y2I1v5VS3W7Sin7oYWinLmckcXzS3ewavsx7m9fjzfvb0Ull9K/\nA96xlEvMiYrns02HuZCeSefg6oztHkyPW2ripLdsVapC00JRDhlj+GjtQd7/4QChAdWYNjSEGmV0\n97tzlzNYvDmRTyPjOJ56mUY1vRnTLYg+beuVeK8spZR90EJRjn214xh/+WI7fpUr8emIDtxSq3KZ\nbTsjK5uvdxxn+rpY9hw/h693JUaEBfBwpwCqebmVWQ6lVOnTQlHObUtMYcy8aC5dyWLSQ+24rUnN\nMt2+MYaoQ2eYERHLz/uT8HB15oFQf0Z1DSKghleZZlFKlQ4tFBXA8dRLjJoTzb4T53i5V3NGhgda\ncsHc/hPnmRkRy4ptR8nMNtzdojZjugfTvkG1Ms+ilCo5WigqiLQrmTy9eBvf7j7JQ50a8HrvFrha\ndHvUk+cuMzcqngUbEzh3OZPQgGqM6R7Mnc1q4awnvpUqd7RQVCDZ2YZ/f7efKT8fIqxhDaY83J6q\nntadL7iYnskX0YnMWh/HkbOXCPL1YlTXIAaE+OuJb6XKES0UFdCymCO8+OVO6lXzYNbwUIL9vC3N\nk5mVzZrdJ5i+LpYdR1Kp7uXG0M4BDO0SgG8Z9dZSShWfFooKakt8Mo/OjyEr2/DJw+0Ja+RrdSSM\nMWyOS2ZGRCw/7D1FJRcn+ofknPhuaHExU0rlTwtFBZaYnMaouVuITbrI3/u05KFODayO9LuDpy4w\na30sy37Nuevenc1qMbZ7MKEB1XTkWqXsjBaKCu785Qwe/3wrP+9PYmR4IC/3am5XJ5STzqczf0M8\n8zYmkJKWQdv6VRnbPZieLWrbVU6lHJkWCgeQmZXNG6v38WlkHD2a+DFpcDsqu9vXaLCXrmSxNCaR\nmevjSDiTRv3qHozuGswDof54urlYHU8ph6aFwoEs3JTAayt3E+znxazhHexyJNisbMP3e04yIyKW\nmISz+Hi4MqRzA4Z3CSyVe3EopW5MC4WDiTp4mnELf8XZSZg2NIQOgdWtjpSvmIRkZqyL49s9J3B1\ncqJvu7qM6RZM4zIcqkQppYXCIcUmXWD03GiOnL3Em/e3on+Iv9WRChR/+iKz1sexJCaRyxnZ3NbE\njzHdg+kSXENPfCtVBrRQOKjUtAzGLYwh6tAZxvVoyHN/amL3w4UnX7zCgo0JzI2K58zFK7SsV4Ux\n3YK5p1Udy65CV8oRaKFwYBlZ2by2ajefbTpMzxa1eH9g23Jx4vhyRhbLtx5lRkQssUkXqVfVg5Hh\ngQzq2ADvSvafX6nyRguFgzPGMDsynn9+vYemtaswa0QodXw8rI5VKNnZhh/3nWJ6RCyb45Kp7O7C\nQ50aMDIsiNo+euJbqZKihUIB8NP+Uzz+2VY83JyZOSyUNvXL161OtyWmMCMilm92HsdJhN5tc058\nN6tTxepoSpV7WijU7w6cPM8jc7aQdD6d9x5sw72t61odqcgSk9P4NDKOxVsSSbuSRbfGvoztHkzX\nRr564lupYtJCof7gzIV0HlsQw5b4szx1Z2OevKNxufyCTU3LYOHmBGZHxpN0Pp2mtSszplsw97Wp\ni5uLnvhWqii0UKjrpGdm8dKXu1j26xHua1OXdwe0LrfDgqdnZrFq2zFmRMRy4OQFalWpxMjwIAZ3\nbICPh31dna6UvdJCofJkjGHqL7G88+0+2vhXZfqwEGpWLr8niI0x/HIgiRkRsUQePIOXmzODOjZg\nZHgg/tXs7wp1peyJFgpVoG93n+CpRduo5unKzOEdaF63/J8c3nU0lZkRsfx3x3EAerWqw9juwbSs\n52NxMqXsU2ELRaEP6oqIs4hsFZGvbK8Xish+EdklIp+KSJ77+yKyRkRSrrbLY/4kEbmQ6/UIEUkS\nkW22x+jCZlSF17NFbZY81gUDDJgaxfd7Tlod6aa1rOfDB4Pase7523gkPJAf953i3knrGTx9Iz/t\nO0V2dvn/UaSUFYpy9u9JYG+u1wuBpkArwAPI7wv9XWBoXjNEJBTIq7/mYmNMW9tjZhEyqiJoWc+H\nlRPCaVzTm7Hzo5n6yyEqwh5mvaoe/LVXc6JevJ2X7mlK3OmLjJyzhZ4frOOLLYmkZ2ZZHVGpcqVQ\nhUJE/IFewO9f2saY1cYG2AzkObCQMWYtcD6PdTqTU0SeL0ZuVUJqVnFn8aNduKdVHd76Zh/PLd1R\nYb5Iq7i7MrZ7Q9Y9fxvvD2yDi7MTzy/bQfhbPzH5p4OkpF2xOqJS5UJh9yg+IOcLPfvaGbZDTkOB\nNUXc9kRglTHmeB7z+ovIDhFZKiL1i7heVUTurs58PLgdT97RmKUxRxg6czPJFyvOl6ibixP92vmz\n+omuLBzdiRZ1q/Dut/vp8uaP/G3Vbg6fSbM6olJ27YaFQkTuBU4ZY2LyWWQKsM4YE1HYjYpIXeAB\nYFIes/8LBBpjWgM/AHPzWcdYEYkWkeikpKTCblrlQ0R4+q5b+GhwO7YfSaHP5PX8dvK6HcFyTUQI\nb+TL3Ec6suapbvRqXYeFmxLo8e+fmLDwV7YePmt1RKXs0g17PYnIm+TsMWQC7kAV4EtjzBAReQ1o\nB9xvjLlubyPXOnoAzxpj7rW97gXMAi7bFmkAxBpjGl3TzhlINsYU2G1Fez2VrK2HzzJmXgzpGVlM\neqgdPZrUtDpSqTl57jJzouJZsDGB85cz6RBYjTHdgrmzWS27H3VXqZtVKt1jc3/h23ojPQLcYYy5\nVNh2+cy/YIzxtj2vc/VwlIj0A14wxnQuaP1aKEresZRLjJobzf4T53j13uYMDwssl1dyF9aF9Ey+\n2JLIrPVxHE25RLCvF6O6BdG/vX+5vShRqRspi0KRCSTwvxPVXxpj/m7ryfSYMWa0rU0EOb2jvIEz\nwChjzLfXrDd3oXgT6E3OHkwyMM4Ys6+gXFooSsfF9EyeWryN7/ecZEjnBrx2X4sKf3+IzKxsvtl1\ngunrYtl5NJXqXm4M6xLA0M4B1PCuZHU8pUqUXnCnSkR2tuHtb/cx7ZdYujbyZfJD7fHxrPhDZBhj\n2BSXzIx1sazdd4pKLk4MCPFnVNcggv28rY6nVInQQqFK1JLoRF5avpP61T2ZNbwDQb5eVkcqMwdP\nnWdmRBxf/nqUjOxs7mpWi7HdgwkJqFahD8epik8LhSpxm+OSeXR+NNkGPhnSnrCGvlZHKlOnzl9m\n/oYE5m9MICUtg3YNqjK2WzB/alEbZz3xrcohLRSqVBw+k8aouVuIO32Rf/RtyeCODayOVObSrmSy\nNOYIMyPiOJycRoPqnozuFsSAEP9ycctZpa7SQqFKzbnLGTz+2VZ+OZDEqK5BvHRPM4f8RZ2Vbfhu\n9wmmR8Sy9XAKVT1dGdo5gGFdAvGrrCe+lf3TQqFKVWZWNv/8ei9zouK5vWlNPhzUlsruFf8kd35i\nEpKZvi6W7/acxNXJiYc6NeDFe5pSyUW71ir7pYVClYkFGxN4bdVuGvl5M3N4KPWrO/Y9IOJOX2T6\nukN8vjmR0IBqTBsaot1qld0q8WHGlcrLkM4BzB3ZkeOpl+g7OZLo+GSrI1kqyNeLN+9vzccPtWPn\n0VT6TI5k/4mKNRSKcjxaKNRN69rYl+UTwqns7sJDMzbx5a9HrI5kuXtb1+WLR7uQnplN/0+i+Gnf\nKasjKVVsWihUiWjo582KCeGEBFTjmS+2886afQ5/o6A29auyamI4Dap7MmruFmatj6sQ9/tQjkcL\nhSoxVT3dmDeqI4M71mfKz4cYv/BX0q5kWh3LUnV8PFg6rgt3NqvFP77aw19X7CIjK9/xM5WyS1oo\nVIlydXbijX6teOXe5ny35wQPTtvAidTLN25YgXm6uTB1SAjjejTks02HGTF7M6lpGVbHUqrQtFCo\nEicijOoaxMzhocSfTqP3x+vZcSTF6liWcnISXri7Kf9+oA2b45LpNyWSuNMXrY6lVKFooVCl5vam\ntVg2LgxXZycenLaB1TvzupmhYxkQ4s9nYzqTcimDvpMjiTp42upISt2QFgpVqprUrszKieG0qOvD\n+IW/Mmntbw5/QrdDYHVWjA+nZuVKDPt0M59tOmx1JKUKpIVClTpf70osHN2Jfu3q8d73B3hq8TYu\nZ2RZHctSDWp4smx8GOGNfHlp+U7+/t89ZDl4LzFlv7RQqDLh7urMfx5sw3M9m7By2zEGTd/IqfOO\nfZK7irsrs4aHMiIskE8j4xgzL5rzl/Ukt7I/WihUmRERJtzWiKlD2rP/xHn6fhzJ3uPnrI5lKRdn\nJ/7WuwX/7NuSXw4kMeCTDSQmp1kdS6k/0EKhytzdLeuw5LEuZBvo/0kUP+w5aXUky107FEpMgmMP\nhaLsixYKZYmW9XxYOTGcRjW9GTM/munrDjn8Se7cQ6EMnr6J5Vt1KBRlH7RQKMvUquLO4rFduKdl\nHd5YvY8Xlu3gSqZjX7Xc0M+b5ePDaR9QlacXb+fdb3UoFGU9LRTKUh5uzkwa3I4nbm/EF9FHGDJr\nE8kXr1gdy1LVvNyY90gnBnWoz+SfDjHhMx0KRVlLC4WynJOT8MyfmvDhoLZsS0yh7+RIDp5y7KG5\n3VycePP+VrzcqxlrdutQKMpaWiiU3ejTth6LxnYm7UoW/aZEse5AktWRLCUijO4WzMxhocQlXaTP\n5PXsPJJqdSzlgLRQKLvSvkE1Vk4Mp15VD0bO2cLcqHirI1nujma1WDouDBcnJx6YFsU3OhSKKmNa\nKJTdqVfVg6XjwritiR+vrdrNKyt2kengQ3M3q1OFFRPCaV6nCuMW/srHP+pQKKrsaKFQdsm7kgvT\nhobyaPdg5m9MYOScLaRecuyrlv0qV+KzMZ3p07Yu//7uAM98sd3hh0JRZUMLhbJbzk7Ci/c0453+\nrdkYe4b7p0QS7+BDc7u7OvPBwLb85a5bWL71KA/P3MTpC+lWx1IVXKELhYg4i8hWEfnK9nqhiOwX\nkV0i8qmIuObTbo2IpFxtl8f8SSJyIdfrSiKyWEQOisgmEQks2j9JVTQPdqjP/FGdOHPxCn2nRLIx\n9ozVkSwlIjx+R2MmP9Se3cdS6fNxJPtPOHYvMVW6irJH8SSwN9frhUBToBXgAYzOp927wNC8ZohI\nKFD1msmjgLPGmEbA+8DbRcioKqjOwTVYOSGcGl5uDJm5icVbdGjuXq3r8MWjXcjIyqb/J1H8tO+U\n1ZFUBVWoQiEi/kAvYObVacaY1cYG2Az459XWGLMWuO7njog4k1NEnr9mVh9gru35UuAOEZHC5FQV\nW0ANL74cH06XhjV4YdlO/vW1Ds3d2r8qKyeGE1DDk1FztzAzIlZPcqsSV9g9ig/I+UK/ruuJ7ZDT\nUGBNEbc9EVhljLm2r189IBHAGJMJpAI18tjuWBGJFpHopCTH7m/vSHw8XJk9ogPDuwQwIyKOsfOi\nuZDu2Fct1/HxYMljXfhT89r88+u9vLR8JxkO3ktMlawbFgoRuRc4ZYyJyWeRKcA6Y0xEYTcqInWB\nB4BJec3OY9p1P5GMMdONMaHGmFA/P7/CblpVAC7OTrzepyX/6NOCnw8kMeCTKI6cdeyhuT3dXJjy\ncHvG92jI55sTGTZrMylpjj0Uiio5hdmjCAd6i0g8sAi4XUQWAIjIa4Af8EwRt9sOaAQctK3XU0QO\n2uYdAerb1u8C+AA65rK6ztAugcwZ2YGjKVeH5j5rdSRLOTkJz9/dlPceaENMwln6TYkiNunCjRsq\ndQM3LBTGmBeNMf7GmEBgEPCjMWaIiIwGegKDjTFF2s81xnxtjKltjAm0rTfNdvIaYBUw3PZ8gG17\netBV5albYz+Wjw/Hq5ILg2dsZMXWo1ZHslz/EH8WjulE6qUM+k6OJPLgaasjqXLuZq6jmArUAjaI\nyDYReRVyejKJyO8nvUUkAlhCzknpIyLS8wbrnQXUsO1hPAP8301kVA6gUU1vVowPp139qjy1eBvv\nfbff4Yfm7hBYnZUTwqnt487wTzfz2SbtJaaKTyrCj/XQ0FATHR1tdQxlsSuZ2byyYheLoxO5p1Vt\n3nugLR5uzlbHstT5yxk8/vlWft6fxMjwQF7u1RxnJ+1EqHKISIwxJvRGy+mV2arCcHNx4q3+OUNz\nf7NLh+YGqOzuysxhoTwSHsTsyHhGzd3C+cuOPRSKKjotFKpCyT00d2zSBR2am5xeYq/e15x/9WtJ\nxG+n6f9JFInJjt1LTBWNFgpVIenQ3Nd7uFMA8x7pyInUy/SZHEl0vHYmVIWjhUJVWDo09/XCG/my\nfEI4Ph6uPDRjE1/+esTqSKoc0EKhKrSrQ3P31aG5f9fQz5vl48MICajGM19s5501+xy+l5gqmBYK\nVeG5uzrz/sC2PPunnKG5H5qxkaTzjj00d1VPN+aN6sjgjvWZ8vMhxi/8lbQrjj0UisqfFgrlEESE\nibc3ZsrD7dlz/Bx9J0ey78Q5q2NZytXZiTf6teKVe5vz3R7tJabyp4VCOZR7WuUMzZ2ZnU3/KVGs\n3XvS6kiWEhFGdQ1i5vBQ4pIu0vvj9ew4kmJ1LGVntFAoh9PavyorJ3Ql2M+b0fOidWhu4PamtVg2\nPgxXZycenLaB1dpLTOWihUI5pNo+7nzxaBfubpEzNPeLX+7kSqZjD83dtHYVVk7M6SU2XnuJqVy0\nUCiH5eHmzOSH2jPxtkYs2pLIsE83cfaiYw/N7ev9x15iTy/e5vC9xJQWCuXgnJyEZ3s24f2Bbfg1\nIYV+UyI5eMqxh+a+2kvsuZ5NWLHtmPYSU1oolALo186fz8d25kJ6Jv2mRBLxm2PfNVFEmHBbIz7R\nXmIKLRRK/S4koBorJoRT18eDEbO3MH9DvNWRLPfnVnVY8miY9hJzcFoolMrFv5ony8aH0eMWP15Z\nuZvXVu4i08HvP93K34eVE7oS5OelvcQclBYKpa7hXcmF6cNCGdMtiLkbEnhkbjTnHHxo7qu9xHo2\nz+kl9tJy7SXmSLRQKJUHZyfhr72a83b/VkQdPM39U6JIOHPR6liW8nRzYcrD7ZlwW0M+35zTSywl\nzbF7iTkKLRRKFWBghwbMH9WJ0xfS6Ts5kk2xZ6yOZCknJ+G5nk1/7yXWd3Ikh5Icu5eYI9BCodQN\ndGlYgxXjw6nm5caQWZv4IjrR6kiWy+kl1onzlzPpNzmS9b+dtjqSKkVaKJQqhEBfL5aPC6dTUA2e\nX7qDN1bvJcvBh+YOCajOignh1PZxZ/jszSzYmGB1JFVKtFAoVUg+nq7MHtmBoZ0DmL4ulkfnx3Ax\n3bGH5q5f3ZNl48Lo3tiXl1fs4m+rdjt8L7GKSAuFUkXg6uzEP/q25PXeLfhx30n6fxLF0ZRLVsey\nVGV3V2YO78Aj4UHMiYpn9DztJVbRaKFQqhiGhwUyZ2RHjp69RJ+PI/n18FmrI1nK2Ul49b7mvNGv\nFet/O03/KVEkJqdZHUuVEC0UShVT91v8WD4hDE83ZwZN38jKbUetjmS5hzo1YN4jHTl1Pp0+kyPZ\nEp9sdSRVArRQKHUTGtWszIoJ4bStX5UnF23jP98fcPj7T4c18mX5+DCqerjy8IxNLIs5YnUkdZO0\nUCh1k6p7ubFgVCceCPHno7W/8fiirVy64thDcwf7ebN8fDihgdX4y5LtvL1mn8MX0PKs0IVCRJxF\nZKuIfGV7vVBE9ovILhH5VERc82m3RkRSrrbLNX2WiGwXkR0islREvG3TR4hIkohssz1G38w/UKmy\n4ObixDsDWvPin5uyeudxBk7fwMlzjn3/aR9PV+Y+0pHBHRvwyc+HGLcwhrQrjt1LrLwqyh7Fk8De\nXK8XAk2BVoAHkN8X+rvA0DxMulr/AAAST0lEQVSmP22MaWOMaQ0cBibmmrfYGNPW9phZhIxKWUZE\nePTWhkwfGsrBUxfo83Eku46mWh3LUq7OTrzRryWv3Nuc7/ec5IGpGzie6ti9xMqjQhUKEfEHegG/\nf2kbY1YbG2Az4J9XW2PMWuB8HtPP2dYt5BQa3S9VFcJdzWux9LEwnAQemLqBNbsc+/7TIsKorkHM\nGt6BhDNp9Pk4ku2JKVbHUkVQ2D2KD4DngeuupLEdchoKrCnqxkVkNnCCnD2TSblm9c91SKp+Uder\nlNWa163CionhNK1TmccW/Mrknw46/NDctzWtybJxYbi5OPHgtA18vcOxC2h5csNCISL3AqeMMTH5\nLDIFWGeMiSjqxo0xI4G65BzSGmib/F8g0HZI6gdgbj65xopItIhEJyU59t3IlH2qWdmdz8d0pneb\nurz77X7+smQ76ZmOfZK7Se2cXmIt6/kw4bNf+Wjtbw5fQMuDwuxRhAO9RSQeWATcLiILAETkNcAP\neKa4AYwxWcBioL/t9RljzNUb9M4AQvJpN90YE2qMCfXz8yvu5pUqVe6uznw4qC3P3HULv5284PDj\nQwH4eldi4ehO9GtXj/98f4CnFm/jcoZjF1B7d8NCYYx50Rjjb4wJBAYBPxpjhth6I/UEBhtjijS4\ni+RodPU5cB+wz/a6Tq5Fe/PHE+hKlTsiwhN3NGbpuC54urlYHccuuLs6858H2/Bczyas3HaMwTM2\nknQ+/cYNlSVu5jqKqUAtYIOtG+urACISKiK/n/QWkQhgCXCHiBwRkZ6AAHNFZCewE6gD/N3W5AkR\n2S0i24EngBE3kVEpu1HJxdnqCHZFRJhwWyM+ebg9e4+fo+/kSPYeP2d1LJUHqQjHB0NDQ010dLTV\nMZRSxbTraCqj5m7hwuVMPhzUjjub17I6kkMQkRhjTOiNltMrs5VSlmtZz4dVE7sS7OfNmPnRzFgX\nqye57YgWCqWUXahVxZ0vHu3Cn1vW5l+r9/J/y3ZyJVPvbWEPtFAopeyGh5szHw9uz8TbGrE4OpGh\nszZx9uIVq2M5PC0USim74uQkPNuzCe8PbMPWwyn0nRLJwVMXrI7l0LRQKKXsUr92/nw+tjMX0zPp\nNyWS9b+dtjqSw9JCoZSyWyEB1VgxIZy6Ph4Mn72Z+RsTrI7kkLRQKKXsmn81T5aND+PWW/x4ZcUu\n/rZqN5lZepK7LGmhUErZPe9KLswYFsrorkHMiYrnkbnRnLucYXUsh6GFQilVLjg7CS/f25w3729F\n1MHT9J8SxeEzaVbHcghaKJRS5crgjg2YN6ojp86n03dKJJvjkq2OVOFpoVBKlTthDX1ZMSGcqh6u\nPDxzI0tjjlgdqULTQqGUKpeCfL1YPj6cjkHVeXbJdt76Zh/ZOox7qdBCoZQqt3w8XZkzsiMPd2rA\n1F8O8diCGC6mZ1odq8LRQqGUKtdcnZ34Z9+WvHZfc37Ye5IHpm7gWMolq2NVKFoolFLlnogwMjyI\nWSM6cDg5jT6TI9memGJ1rApDC4VSqsK4rUlNvhwfRiUXJx6ctoGvdhyzOlKFoIVCKVWh3FKrMisn\nhNOqng8TP9vKhz/8pve2uElaKJRSFU4N70osHNOJ+9vX4/0fDvDkom1czsiyOla5pXd6V0pVSJVc\nnHnvgTY0qunNO2v2czg5jenDQqhZ2d3qaOWO7lEopSosEWF8j0ZMHRLC/hPn6ftxJHuOnbM6Vrmj\nhUIpVeHd3bI2Sx7rQraBAVOj+GHPSasjlStaKJRSDqFlPR9WTgynUU1vxsyPZvq6Q3qSu5C0UCil\nHEatKu4sHtuFe1rW4Y3V+3hh2Q6uZOq9LW5ET2YrpRyKh5szkwa3o6GfFx/9eJCEM2lMHRJCNS83\nq6PZLd2jUEo5HCcn4Zk/NeHDQW3ZmphC3ymRHDx1wepYdksLhVLKYfVpW4/Px3TmYnom/aZEsu5A\nktWR7JIWCqWUQwsJqMaKCeHUq+rByDlbmL8h3upIdqfQhUJEnEVkq4h8ZXu9UET2i8guEflURFzz\nabdGRFKutss1fZaIbBeRHSKyVES8bdMrichiETkoIptEJLD4/zyllLox/2qeLB0XRo9b/Hhl5W5e\nW7mLzCw9yX1VUfYongT25nq9EGgKtAI8gNH5tHsXGJrH9KeNMW2MMa2Bw8BE2/RRwFljTCPgfeDt\nImRUSqli8a7kwvRhoYzpFsTcDQk8Mjeac5czrI5lFwpVKETEH+gFzLw6zRiz2tgAmwH/vNoaY9YC\n5/OYfs62biGn0Fzt0NwHmGt7vhS4w7aMUkqVKmcn4a+9mvN2/1ZEHTzN/VOiSDhz0epYlivsHsUH\nwPPAdftitkNOQ4E1Rd24iMwGTpCzZzLJNrkekAhgjMkEUoEaebQdKyLRIhKdlKQnoJRSJWdghwbM\nH9WJ0xfS6Ts5kk2xZ6yOZKkbFgoRuRc4ZYyJyWeRKcA6Y0xEUTdujBkJ1CXnkNbAq5vMa9E82k43\nxoQaY0L9/PyKummllCpQl4Y1WDE+nGpebgyZtYkvohOtjmSZwuxRhAO9RSQeWATcLiILAETkNcAP\neKa4AYwxWcBioL9t0hGgvm39LoAPkFzc9SulVHEF+nqxfFw4nYJq8PzSHbz5zV6ysx1v2I8bFgpj\nzIvGGH9jTCAwCPjRGDNEREYDPYHBxpgidQ+QHI2uPgfuA/bZZq8ChtueD7Btz/H+zyil7IKPpyuz\nR3ZgSOcGTPsllkcXxHAxPdPqWGXqZq6jmArUAjaIyDYReRVAREJF5PeT3iISASwh56T0ERHpSc7h\npbkishPYCdQB/m5rMguoISIHydlT+b+byKiUUjfN1dmJf/Rpyd/ua87avScZMHUDx1IuWR2rzEhF\n+LEeGhpqoqOjrY6hlHIAP+8/xeOfbcXdzZkZw0JpW7+q1ZGKTURijDGhN1pOr8xWSqki6NGkJl+O\nD8Pd1YmB0zawavsxqyOVOi0USilVRI1rVWbF+HBa+/vwxOdbef/7AxX63hZaKJRSqhhqeFdiwehO\n9G/vz4drf+OJRdu4nJFldaxSofejUEqpYqrk4sy/H2hNo5revPPtPhKT05g+LISald2tjlaidI9C\nKaVugogwrkdDpg4JYf+J8/T9OJI9x85ZHatEaaFQSqkS0LNFbZY81gUDDJgaxfd7TlodqcRooVBK\nqRLSsp4PKyeE07imN2PnRzP1l0MV4iS3FgqllCpBNau4s/jRLtzTqg5vfbOP55fu4Epm+b63hZ7M\nVkqpEubu6sykQe1o6OfNR2t/IyE5jalDQqju5WZ1tGLRPQqllCoFTk7CM3fdwoeD2rItMYW+kyM5\neOq6W/OUC1oolFKqFPVpW49FYzuTdiWLflOiWHeg/N0/RwuFUkqVsvYNqrFyYjj1qnowcs4W5kbF\nWx2pSLRQKKVUGahX1YOl48K4rYkfr63azSsrdpGZVT5OcmuhUEqpMuJdyYVpQ0MZ2z2Y+RsTGDln\nC6mXMqyOdUNaKJRSqgw5Owkv3dOMt/u3YsOhM9w/JZKEMxetjlUgLRRKKWWBgR0asGB0J85cvELf\nyZFsij1jdaR8aaFQSimLdA6uwYrx4VT3cmPIrE18EZ1odaQ8aaFQSikLBfp68eX4cDoH1+D5pTt4\nY/VesrLta9gPLRRKKWUxHw9XZo/owNDOAUxfF8uj82O4mJ5pdazfaaFQSik74OLsxD/6tuT13i34\ncd9JBkzdwNGUS1bHArRQKKWUXRkeFsjskR05kpxGn48j2Xr4rNWRtFAopZS9ufUWP74cH4anmzMD\np29k1fZjlubRQqGUUnaoca3KrJgQTlv/qjzx+Vb+8/0By+5toYVCKaXsVHUvN+aP7siAEH8+Wvsb\nEz/fyuWMrDLPofejUEopO1bJxZl3B7SmcU1v3lqzjyPJacwYFkrNKu5llkH3KJRSys6JCI/e2pBp\nQ0I4cPICfSZHsvtYapltv9CFQkScRWSriHxle71QRPaLyC4R+VREXPNpt0ZEUq62yzU9z/Yi0kNE\nUkVkm+3x6s38A5VSqqL4U4vaLB3XBYAHpm7gu90nymS7RdmjeBLYm+v1QqAp0ArwAEbn0+5dYGge\n0wtqH2GMaWt7/L0IGZVSqkJrUdeHlRPCaVyrMo8uiOHT9XGlvs1CFQoR8Qd6ATOvTjPGrDY2wGbA\nP6+2xpi1wHX3/ytse6WUUn9Us4o7i8d2pnebugT5eZX69gq7R/EB8Dxw3V02bIeMhgJrihMgn/Zd\nRGS7iHwjIi2Ks16llKrI3F2d+XBQO25rUrPUt3XDQiEi9wKnjDEx+SwyBVhnjIkoZoZr2/8KBBhj\n2gCTgBX55BorItEiEp2UVP7uQauUUuVFYfYowoHeIhIPLAJuF5EFACLyGuAHPFOcjefV3hhzzhhz\nwfZ8NeAqIr7XtjXGTDfGhBpjQv38/IqzeaWUUoVww0JhjHnRGONvjAkEBgE/GmOGiMhooCcw2BhT\n5Bu/5tdeRGqLiNied7RltN87eiilVAV3M9dRTAVqARtyd2MVkVAR+f2kt4hEAEuAO0TkiIj0LKg9\nMADYJSLbgY+AQcaq69aVUkohFeE7ODQ01ERHR1sdQymlyhURiTHGhN5oOb0yWymlVIG0UCillCqQ\nFgqllFIFqhDnKEQkCUiwvfQFTlsYp6g0b+nSvKVL85au0s4bYIy54fUFFaJQ5CYi0YU5OWMvNG/p\n0rylS/OWLnvJq4eelFJKFUgLhVJKqQJVxEIx3eoARaR5S5fmLV2at3TZRd4Kd45CKaVUyaqIexRK\nKaVKkN0WChFpKyIbbeNARdsGCEREnst1m9RdIpIlItXzaD9RRA6KiMlr9FkR6WBrOyDXtAYi8p2I\n7BWRPSISaM95bdOriMhREfm4sFmtyGvb3gYR2S0iO0RkoD3ntU0bLiK/2R7D7SGviPSx/f2urrdr\nrnnv2P6+e0XkI5GcwTXtOG+xP29WZbbNL/Jnrqyz3uzn7TrGGLt8AN8Bf7Y9vwf4OY9l7iNnNNu8\n2rcDAoF4wPeaec7Aj8BqYECu6T8Dd9meewOe9pzXNu9D4DPgY3v++wK3AI1tz+sCx4Gqdpy3OhBr\n+2812/NqVue1vS+vHjJuDeyzPQ8DIm3/FmdgA9DDXvPe7OfNqszF/cxZ8H64qc/btQ8X7JcBqtie\n+wDH8lhmMPB5no2N2QqQz4+qx4FlQIerE0SkOeBijPne1v6CPee1LRtCzgi8a4Ci9rUu07zGmAO5\nnh8TkVPk3IskxR7zkjME/vfGmGRbu++Bu/Nbf1nlveZ96WXbztXtuQNugACuwMlCZi3zvCXweSvz\nzLZli/uZK9OsJfB5uy6AXT6AZsBhIBE4Ss4VhLnnewLJQPUbrCeeP1bgesAv5PzqmsP/fkH2Bb4C\nvgS2Au8Cznac14mcX2T1gREUfY+iTPNe06YjsBdwste8wLPAy7mWewV41uq8tmn9gH229l1yTf83\nOV8EqcC/7OH9kF9ebvLzZlHmYn/mrHg/5Jpf5M/btQ9L9yhE5Aegdh6z/grcATxtjFkmIg8Cs4A7\ncy1zHxBpbL/4iuAD4AVjTNY11dkF6EbOLt5hYDE5b4ZZdpp3PLDaGJOY36FoO8t7NVMdYD4w3Fxz\nwys7y5vXH/UPXQQtyosxZjmwXES6A/8A7hSRRuR8GfnbFvteRLobY9bZY14K8Xmzw8wFfubsLOvV\nTPl+3oq6Ebt8kPOr6OqxNwHOXTN/OfBQIdYTzx9/QcbZpsUDF4BT5Py66Uyu44bAUGCyHeddSM4H\nLJ6csWDOAW/Za17bvCrk3BP9gXLwfhgMTMu13DRy7sZoad485seRMx7Qc8Aruaa/Cjxvx3lv6vNm\nUeZif+bKOqvtebE/b9et92ZXUFoPcnaVetie3wHE5JrnQ85ultfN/GH546EGZ2A74Gd7PRuYYK95\nr5k+gqIfeirrv68bsBZ4qpy8H6rbPnTVbI84bnBYoCzyAo1yfeG0J+cwhgADgR/I+aXuavtb32fH\neW/q82ZF5mvaFOkzZ8Hf96Y+b9dttyRWUhoPoCsQY3szbQJCrvmftCiPNquBurbnTwBHgExyThzN\nzGP5Ofyx19NdwA5gp22emz3nLe6b1oq8wBAgA9iW69HWXvPaXj8CHLQ9RtrD3xd4Adht+/ttALra\npjuTs9ezF9gD/Mee897s582qzMX9zFnwfripz9u1D70yWymlVIHs9oI7pZRS9kELhVJKqQJpoVBK\nKVUgLRRKKaUKpIVCKaVUgbRQKKWUKpAWCqWUUgXSQqGUUqpA/w/mAF4NbDkg/AAAAABJRU5ErkJg\ngg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x69876122e8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "test_oco_box = area[0]\n",
    "#1432\n",
    "p1 = (test_oco_box[0][0], test_oco_box[1][0])\n",
    "p2 = (test_oco_box[0][3], test_oco_box[1][3])\n",
    "p3 = (test_oco_box[0][2], test_oco_box[1][2])\n",
    "p4 = (test_oco_box[0][1], test_oco_box[1][1])\n",
    "\n",
    "plt.plot(area[0][0], area[0][1])\n",
    "print(area[0])\n",
    "test_polygon = Polygon([p1, p2, p3, p4])\n",
    "    \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-87.1454315186 42.1253089905\n",
      "-87.1454315186 42.1293089905\n",
      "-87.1454315186 42.1333089905\n",
      "-87.1454315186 42.1373089905\n",
      "-87.1454315186 42.1413089905\n",
      "-87.1454315186 42.1453089905\n",
      "-87.1454315186 42.1493089905\n",
      "-87.1394315186 42.1253089905\n",
      "-87.1394315186 42.1293089905\n",
      "-87.1394315186 42.1333089905\n",
      "-87.1394315186 42.1373089905\n",
      "-87.1394315186 42.1413089905\n",
      "-87.1394315186 42.1453089905\n",
      "-87.1394315186 42.1493089905\n",
      "-87.1334315186 42.1253089905\n",
      "-87.1334315186 42.1293089905\n",
      "-87.1334315186 42.1333089905\n",
      "-87.1334315186 42.1373089905\n",
      "-87.1334315186 42.1413089905\n",
      "-87.1334315186 42.1453089905\n",
      "-87.1334315186 42.1493089905\n"
     ]
    }
   ],
   "source": [
    "#getting fpar\n",
    "\n",
    "from my_functions import coords_to_ind\n",
    "\n",
    "\n",
    "import numpy as np\n",
    "for i in np.arange(min(test_oco_box[0]), max(test_oco_box[0]), 0.006):\n",
    "    for j in np.arange(min(test_oco_box[1]), max(test_oco_box[1]), 0.004):\n",
    "        \n",
    "        print(i,j)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "gnome_test_data = Dataset('08\\\\ret_f_nr5_nsvd12_v26_waves734_nolog.20160801_v27_all.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-162.73866272, -161.11044312, -163.13847351, -161.51272583],\n",
       "       [-162.73866272, -161.11044312, -163.41612244, -161.79342651],\n",
       "       [-163.3377533 , -161.71464539, -164.04302979, -162.42326355],\n",
       "       ..., \n",
       "       [ 145.1022644 ,  145.26411438,  142.3263855 ,  142.50340271],\n",
       "       [ 142.40600586,  142.58261108,  139.69844055,  139.89016724],\n",
       "       [ 139.76766968,  139.95896912,  137.08872986,  137.29917908]], dtype=float32)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cur_data.variables['Longitude_corners'][:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 75.94194031,  76.0980835 ,  76.21187592,  76.36787415],\n",
       "       [ 75.94194031,  76.0980835 ,  76.3137207 ,  76.47038269],\n",
       "       [ 76.26849365,  76.42510986,  76.63511658,  76.79251862],\n",
       "       ..., \n",
       "       [-58.29808807, -57.87671661, -57.99367142, -57.6099968 ],\n",
       "       [-58.00455093, -57.61994553, -57.65233612, -57.29317474],\n",
       "       [-57.66324997, -57.30360031, -57.25914383, -56.91117477]], dtype=float32)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cur_data.variables['Latitude_corners'][:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0xf27b926dd8>]"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3XtwnPd93/v3D/c7sLsASIIgsAuA\npMSLxAsASiQl62KrsnxNznEixYnt+CRy0kp12/G4dd06nU49deszp06azLQaVzn19MRxpid23CiJ\nJrajJBIlAbyKpCiaBAGQ4B3YxY0gbrvf/vE8fAhRIAmSAHaB/bxmOMI+u8/uF6sH3+fy+/6+jzMz\nREQke+SkOwAREVlcSvwiIllGiV9EJMso8YuIZBklfhGRLKPELyKSZZT4RUSyjBK/iEiWUeIXEcky\neekOYDbV1dUWjUbTHYaIyJKxb9++fjOrmctrMzLxR6NR9u7dm+4wRESWDOdc71xfq0s9IiJZRolf\nRCTLKPGLiGQZJX4RkSyjxC8ikmWU+EVEsowSv4hIlsnIOn4RyV6T0ykOnx2isyfO5tWV7GqpTndI\ny44Sv4ik1djkNAdOD/J2d5zO7jgHziQYn0oB8NuPNSvxLwAlfhFZVINjk+ztSdDRE6ejO86Rs0NM\np4wcBxvqKniuvYEdsTCt0TDVZYXpDndZUuIXkQV1cXicjm4vyXf2xHnvwggABbk5PLimkucfbaI9\nFmZ7Y4jyovw0R5sdlPhFZN6YGb0DY8HRfGdPnN6BMQBKCnLZ3hjiY5tX0R4L8+CaKoryc9MccXZS\n4heRu5ZKGccvjtDZEw+u0V8amQAgVJJPWzTMrz3USHsszIZVFeTlqpAwEyjxi8icTSX9ihv/0s3e\n3gRDV6cAWFVZxMPNEdqiYXbEwjTXlJGT49IcscxGiV9EburqZJIDp68PxB44PcjVqSQATTWlfHTT\nStpjYdqiYepDxTinRL8UKPGLSGBobIq9vfEg0R/u8ypunIMNqyr45bY1QcVNTbkqbpYqJX6RLHZp\neJyOHu/a/NvdcY5fHMEM8nMdD9ZX8ZszKm4qVHGzbCjxi2QJM+N0fOx9pZU9N1TcPONX3GxRxc2y\npsQvskylUsbPL40ER/OdPXEuDnsVN1V+xc1nd/gVN3UV5KviJmso8YssE1PJFEf8HjfeEf31ipuV\nFUXsiEVoi3kVNy2quMlqSvwiS9TVySQHziTo7E7Q0TPA/t7rFTex6lKe3rgySPSquJGZlPhFloih\nq1Ps643T0Z2go3uAw2eHmEp6FTf3r/QqbtqiYdpiIWrLi9IdrmQwJX6RDHVpZJzO7kQwK/a9C8NB\nxc3m1ZX8X7ub2BELs60xRGWxKm5k7pT4RTKAmdGXuBq0PejoidPdfwWA4vxctjVW8U+eXEdbLMTW\nNSGKC1RxI3dPiV8kDVIp48Sl0aCGvqM7zoXhcQAqi72Km+favUs3m1ZXquJG5pUSv8gimE6mOHpu\n2Kuh7/FKKwfHvIqbFRWFtMcitEdDtMcirK1VxY0sLCV+kQUwPpXk4JnBYKLUvt4EY5NexU00UsJT\nG1b4zcwirAmr4kYWlxK/yDwYHp9in39Xqc7uOIf6BoOKm/UryvnM9nraYmHao2FqK1RxI+mlxC9y\nFy6PTAQTpTr8ipuUQV6OY3N9JV/cHaM9Gqa1MUxliSpuJLMo8YvcxrWKm2uXbTq645zyK26K8nPY\n1hDiHz+5lvZomC0NVZQU6M9KMpu2UJEbmBknL40G/W06uuOcH/IqbiqK8miLhr3JUrEwm+oqKchT\nxY0sLbdN/M659cAPZixqAr5hZt9xzr0IvABMA6+Y2Vdv8h65wF7grJl9/N7DFpk/08kU754ffl/X\nyoRfcVNbXhi0PWiLhlm/olwVN7Lk3Tbxm9lxYAsECfws8EPn3OPAp4AHzGzCOVd7i7f5MnAMqLj3\nkEXuzfhUkkN+xU1HT5z9vQmu+BU3jZESnrx/Be3+QGxjpEQVN7Ls3OmlnieBLjPrdc59G/iWmU0A\nmNml2VZwztUDHwO+CfyzewlW5G6MjE+xrzcRHM0fOjPEZDIFwH0ry/nFbfVeoo+FWaGKG8kCd5r4\nnwW+7/+8DnjEOfdNYBz4ipl1zrLOd4CvAuV3HaXIHegfnWBvz/Ue9O+eu15xs2l1JV/YFfUqbqIh\nqkoK0h2uyKKbc+J3zhUAnwS+NmPdEPAQ0Ab8iXOuycxsxjofBy6Z2T7n3GO3ef/ngecBGhoa7uR3\nkCzXlxh7X2ll1+XrFTdb14R48Ym1tMfCbFXFjQhwZ0f8HwX2m9lF/3Ef8Kd+ou9wzqWAauDyjHV2\nAZ90zj0DFAEVzrn/YWa/euObm9lLwEsAra2tduPzIuBV3HRdHr3ezKw7zjm/4qbcr7j5TKvX42bz\nalXciMzmThL/c1y/zAPwI+AJ4DXn3DqgAOifuYKZfQ3/DME/4v/KbElf5GamkymOnR/h7e4BOnu8\nu0rFr0wCUFNeSHs0zJeuVdysLCdXFTcitzWnxO+cKwE+AnxpxuKXgZedc0eASeDzZmbOuTrgu2b2\nzLxHK8ve+FSSd/qG6OgeoKMnwb6eeFBx0xAu4fH1tV5pZSxMVBU3IndlTonfzMaAyA3LJoHZLtmc\nAz6Q9M3sNeC1uwlSlq/RiWm/4maAzu4EB88MBhU361eU8wvbVvudK8OsrFTFjch80EiXLKqB0Qk6\ne66XVh49N0TKINevuPn8zkbaYxFaG0OESlVxI7IQlPhlQZ0dvEpn9/XSypOXRgEozMtha0MVLzze\nQnsswtaGKkoLtTmKLAb9pcm88SpurryvtPLs4FUAygvzaI2G+MVtq9kR8+4qVZin2weKpIMSv9y1\nZMo4dkOPmwG/4qa6rJD2WIjfeCRGeyzMfSsrVHEjkiGU+GXOJqavVdx4iX5fb4LRiWkA6kPFfGh9\nTdDMLFZdqoobkQylxC83NToxzX6/x01HT9yruJn2Km7W1pbxqS11tPuJvq6qOM3RishcKfFLIH5l\nMrg+71XcDJNMmVdxU1fB5x5qpD0WpjUaJqyKG5ElS4k/i50bvErntWZm3XFO+BU3BXk5bF1TxT98\nrNnvcROiTBU3IsuG/pqzhJlxqv9K0N+moydOX8KruCnzK24+vdWruNlcr4obkeVMiX+ZulZxM/PS\nTf+oV3ETKS2gPRbmi7u8ipv7V6niRiSbKPEvExPTSQ73DdHhJ/p9PQlG/Iqb1VXFPLq2hjb/ZiNN\nqrgRyWpK/EvUlYlp9p9OBLNiD54ZZMKvuGmpLeMTW+poj3rNzFar4kZEZlDiXyISN1TcHPErbnIc\nbKyr5FcfaqQtGqYtGiJSVpjucEUkgynxZ6jzQ1eDJN/RHefnF69X3Gypr+K3P9RMWyzMtoYqyovy\n0xytiCwlSvwZwMzo7r/W4yZBR88AZ+LXK262N4b41JbVtEXDPFBfSVG+Km5E5O4p8adBMmW8d2HY\nK630k33/6AQA4dIC2qIhvrAzxo5YmPtWlpOXq9sHisj8UeJfBJPTKQ6fHfSO5rsH2NubYGT8esXN\n7paId7ORWIjmmjJV3IjIglLiXwBjk9Ps7x30j+YHOHhmkPEpr+KmuaaUjz9QR3ssRFs0TH2oJM3R\niki2UeKfB4Njk3T2JIL2B0fPDjHtV9xsqKvgV9obaY+FaI2GqVbFjYikmRL/XbgwNE5HTzxof3D8\n4ggABbk5PLimki99qIm2aJjtjSFV3IhIxlHivw0zo3dgLOhv09Ed53R8DIDSgly2NYb4xIOraIuG\neXBNlSpuRCTjKfHfIJUyjl8ceV+ivzziVdyESvJpi4b53MNee+INqypUcSMiS07WJ/7J6RRHznl3\nler0J0wN+xU3qyqL2NkcoT0Wpj0aprmmjBw1MxORJS7rEv/Y5DQHTg8Gtw88cCYRVNw01ZTyzOZV\nwV2l6kPFKq0UkWVn2Sf+obEpOnviQcXNEb/ixjnYsKqCZ9sa2OHfVaqmXBU3IrL8LbvEf3F4/H09\nbo5fHMHMq7h5oL6S5x9toi3mVdxUqOJGRLLQskr8X/x/O/nZe5eCxw83RfinH15HeyzMFlXciIgA\nyyzxP7K2mgtD47x7fhiAI+eGKC/Ko6okn5ryQt2AREQEcGZ26xc4tx74wYxFTcA3zOw7zrkXgReA\naeAVM/vqDesWAX8HFOLtZP6nmf3O7YJqbW21vXv33tEvMtPA6ARvnhrgjZP9/P2J/uDesl6VTjW7\n10bY1VxNbUXRXX+GiEgmcc7tM7PWOb32don/hjfOBc4CO/B2AF8HPmZmE865WjO7dMPrHVBqZqPO\nuXzgdeDLZvbWrT7nXhP/jU4PjPH6yX7eONnPG139DI5NAbBuRRm7WqrZ3VLNjqYIZYXL6gRIRLLI\nnST+O810TwJdZtbrnPs28C0zmwC4Men7ywwY9R/m+//mvqeZJw2REn4l0sCv7GgglTLePT/MGyf7\nef1kP3/09mn+8I0ecnMcW9ZUBTuCLWuqKMjT5CwRWX7u9Ij/ZWC/mf2+c+4g8GfA08A48BUz65xl\nnVxgH9AC/IGZ/fPbfc58H/HfyvhUkv2nE97ZwMkB3ukbJGVQUpDLjljY2xGsrWb9inKND4hIxlqQ\nSz3OuQLgHLDRzC46544APwO+DLThjQM02U3e0DlXBfwQeNHMjszy/PPA8wANDQ3be3t75xTXfBu6\nOsVb/vjA6yf7OXX5CgDVZQXe+EBLNbvWVusG5iKSURYq8X8K+Edm9pT/+K/wLvW85j/uAh4ys8u3\neI/fAa6Y2f99q89azCP+2zk3eJU3Tvazp2uA10/2B317YtWl7GqJsLulmoebqqks0ZwAEUmfhbrG\n/xzw/RmPfwQ8AbzmnFsHFAD9NwRSA0yZ2aBzrhj4MPAf7uAz066uqpjPtK7hM61rMDNOXBrl9RPe\nQPEP95/lf7x1Gudg8+rKYHxge2NIcwZEJGPN6YjfOVcCnMG7lDPkLysAXga2AJN41/h/5pyrA75r\nZs845x4A/juQC+QAf2Jm//Z2n5dJR/y3MpVM8U7fIK+f8C4N7T+dYDplFObl0BYNs9M/I9hYV0mu\nmruJyAJasHLOxbJUEv+NrkxM09EdD0pH37vg3aClsjifnc2R4IygMVKigWIRmVcLWc4pt1BamMfj\n99Xy+H21AFwemWBPl7cTeP1EP3955ALg3WB9V4u3I9jZXK3mcCKyqHTEv0jMjB5/Itkef7B46Ko3\nkey+leVBtVB7NEypJpKJyB3SpZ4lIJkyjp4bCi4LdfYkmJxOkZfj2NYQ8ucPRHigvop83eVLRG5D\niX8JGp9Ksq83EewIDp8dwgzKCvN4qCns9xiqZm1tmcYHROQDdI1/CSrKz2VXSzW7WqoBGByb5E1/\n7sCergF+cszriFFTXuhdFmqpZldLhFWVmkgmIndGR/xLRF9ijD0nB4IzgoErk4B3u8hrO4KHmiJU\nFmsimUg20qWeZS6VMo5fHAnaSnR0xxmbTJLj4IH6qmBHsK2xisI8TSQTyQZK/FlmcjrFwTODwdnA\nwTODJFNGUb43kezajmDDqgpyNJFMZFlS4s9yI+NTvH0q7o8P9PPzi15n7FBJPjtbqtnlN5triJSk\nOVIRmS8a3M1y5UX5fHjDCj68YQXg3YB+T1d/0FrilXfOA7AmXBycDexsriZcWpDOsEVkkeiIP8uY\nGV2Xr/g7gn7ePDXAyPg0ABvrKoLKovZomOICjQ+ILBW61CNzNp1McfjsUDBQvL93kMlkioLcHLY1\nXh8o3ry6kjxNJBPJWEr8cteuTibp7IkHO4Kj54YBKC/M46HmSLAjaK4p1UQykQyia/xy14oLcnl0\nXQ2PrqsBYGB0gjf9O5K9cXKAv373IgArK4qCSWS7WqpZUVGUzrBF5A7oiF/uyOmBMd7o6g+azSXG\nvEZza2vLgrbTO5rClBdpIpnIYtKlHlkUqZTx7vlhb6D45AAd3QOMT6XIzXFsWVPFLv8eBFsbQhTk\naXxAZCEp8UtaTEwn2d87GIwPvNM3SMqgOD+XHU3XJ5KtX1GuiWQi80yJXzLC0NUp3vbHB14/2U/X\n5SsAVJcVsLP5+vhAfUgTyUTulQZ3JSNUFufz1MaVPLVxJQDnh67yxslrA8X9/PjQOQCikZJgfODh\n5ghVJZpIJrKQdMQvaWFmnLw0GvQXeutUnNGJaZyDTXWVwY6gNRqiKF8TyURuR5d6ZMmZSqZ4p28w\naCtx4EyCqaRRkJdDW9S7I9mu5mo2ra4kV+MDIh+gxC9L3pWJaTp64rxxwhsfeO/CCAAVRXne+MBa\n74wgGinRRDIRdI1floHSwjweX1/L4+trAbg8MsGerv7gZjR/dfQCAKurioNB4p3N1dSUF6YzbJEl\nQUf8suSYGb0DY8H4wJ6uAYauehPJ7ltZHowPtMfClBbq2Eaygy71SFZJpoyj54aCiqGOnjiT0yny\nchxbG6qCHcGDa6rIV6M5WaaU+CWrjU8l2debCM4IDp8dwgxKC3J5qCkStJ5et6JM4wOybOgav2S1\novzcILkDDI5N8tapazeqH+Cn710CoKa8MGgrsaulmrqq4nSGLbJodMQvWacvMcaekwO80eWdEfSP\nTgLQVFPKbn+Q+OHmCJXFajQnS4cu9YjMkZlx/OIIr5/wdgJvd8cZm0yS42BzfRW7/YqhbQ2aSCaZ\nbV4Tv3NuPfCDGYuagG+Y2Xeccy8CLwDTwCtm9tUb1l0DfA9YCaSAl8zsd28XlBK/pMvkdIqDZwaD\nthIHzgySTBlF+Tm0RcPBQPGGVRVqNCcZZcGO+J1zucBZYAfeDuDrwMfMbMI5V2tml254/SpglZnt\nd86VA/uAT5vZu7f6HCV+yRQj41N0dMeDgeKfXxwFIFSS7zea83YEDRE1mpP0WsjB3SeBLjPrdc59\nG/iWmU0A3Jj0/WXngfP+zyPOuWPAauCWiV8kU5QX5fPk/St48v4VAFwaHmdP10CwI3jl8HkA6kPF\nQdvpnc0RImWaSCaZ606P+F8G9pvZ7zvnDgJ/BjwNjANfMbPOW6wbBf4O2GRmw7M8/zzwPEBDQ8P2\n3t7eO/g1RBafmXGq/0pwWWhP1wAj49MAbFhVwe613k6gPRampEAFdLKwFuRSj3OuADgHbDSzi865\nI8DPgC8DbXjjAE02yxs658qAvwW+aWZ/ervP0qUeWYqmkymOnBv27j9wop99vQkmkynycx3bGkLe\nGcHaah5YXUmeJpLJPFuoxP8p4B+Z2VP+47/Cu9Tzmv+4C3jIzC7fsF4+8OfAq2b2/8zls5T4ZTm4\nOpmksyfunRF09XP03DBmUF6Yx46mCLtbIuxeW01zjSaSyb1bqGv8zwHfn/H4R8ATwGvOuXVAAdB/\nQyAO+G/AsbkmfZHlorggl0fX1fDouhoA4lcmeXPG+MBPjl0EYEVFYTBIvKulmhUVRekMW7LAnI74\nnXMlwBm8SzlD/rIC4GVgCzCJd43/Z865OuC7ZvaMc2438PfAYbxyToB/aWZ/cavP0xG/ZIPTA2PB\nJLI9XQPEr3gTyVpqy4KdwI6mMBVFmkgmt6cJXCJLTCplHLvgjw+cHKCje4DxqRS5OY4H6yuDHcHW\nhhAFeRofkA9S4hdZ4iamkxw4PRjcqP7QmUFSBsX5ubTHwsGO4L6V5ZpIJoASv8iyM3R1irdPDQQ7\ngq7LVwCIlBaws6Wa3S0RdjZXsyasiWTZSt05RZaZyuJ8ntq4kqc2rgTgwtB4MH/g9ZP9/K9D5wBo\njJQEA8UPN0UIlRakM2zJUDriF1nizIyTl0aDttNvnRpgdGIa52BTXSU7WyLsbqmmLRpWo7llTJd6\nRLLYdDLFob6h4GzgwOkEU0mjIC+H1sZQcEawaXUluRofWDaU+EUkMDY5TUd3PKgYOnbe65hSUZTH\nw82RYKA4Vl2qiWRLmK7xi0igpCCPx9bX8tj6WgD6RyfY0zXAGye8M4JXj3oTyeoqi7yzgbXejWhq\nyzWRbLnSEb9IFjMzTsfHgtnEe7oGGBybAmD9inJ/RxChPRahrFDHiZlMl3pE5K4kU8a754Z5/WQ/\ne7r66eiOMzGdIi/HsbWhip3N3hnBljVV5KvRXEZR4heReTE+lWR/byI4I3jn7BBmUFqQy46mSDBQ\nvG6FGs2lm67xi8i8KMrPZWdLNTtbqgEYGpvizVNe2egbJ/v52Xve/ZeqywrZ1XJ9R1BXVZzOsOU2\ndMQvInft7ODVYCLZGyf76R/1Gs01VZeyy68WergpQmWJGs0tNF3qEZFFZ2YcvzgSnA28dWqAsckk\nOQ42r64Mzga2NYY0kWwBKPGLSNpNTqc41DfI6ye8s4EDZwZJpozCvBzaY2FvoLilmg11FZpINg+U\n+EUk44xOTNPRPcDrJ7wzguMXRwCoKslnZ/P18YGGcIkGiu+CBndFJOOUFebxxH0reOK+FQBcGh5n\nz4w7kv3F4QsA1IeK2dXs3Z94Z3OE6rLCdIa9LOmIX0TSzszo7r8S9Bfa0zXAyPg0APevqmC3XzHU\nHgtTUqDj1dnoUo+ILGnJlHH47FBQLbS3J8FkMkV+rmNrQyjoL/RgfSV5mkgGKPGLyDJzdTLJ3t54\ncFno6LlhzKC8MI8dTZHgjKClNnsnkukav4gsK8UFuTyytoZH1tYAkLgyyZunro8P/OSY12iutrww\nOBvY1VLNyko1mpuNjvhFZMk7Ex973/hA/Io3kayltoxdfsXQQ80RKoqW70QyXeoRkayVShnvXRgJ\ndgQd3XGuTnkTyR5cUxWcEWxtqKIwb/lMJFPiFxHxTUwnOXB6kD3+juBQ3xDJlFGcn0t7LBz0GLp/\nZQU5S3gimRK/iMhNDI9P8fapeHBGcPLSKADh0gJ2zrgj2ZpwSZojvTMa3BURuYmKonw+smEFH9ng\nTSS7MDTOnq7+YKD4z985D0BDuCSYTbyzOUKotCCdYc8rHfGLiPjMjK7Lo7x+wrs/8VunBhidmMY5\n2FhXEewIWhvDFBdk1viALvWIiMyD6WSKQ31DwfjA/tMJppJGQW4O2xtD7F7rXRbavLoy7Y3mlPhF\nRBbA2OQ0Hd3XxgcGOHZ+GICKojwe9stGd7VU01RduugTyeb1Gr9zbj3wgxmLmoBvmNl3nHMvAi8A\n08ArZvbVWdZ/Gfg4cMnMNs0lKBGRTFRSkMdj62t5bH0tAP2jE7zZ5XUb/fsT/bx61JtItqqy6Pr4\nQEuE2vLMmkh2R0f8zrlc4CywA28H8HXgY2Y24ZyrNbNLs6zzKDAKfG+uiV9H/CKy1JgZp+NjwSDx\nnq4BBsemAFi3oizYEexoilBWOP91NQtZ1fMk0GVmvc65bwPfMrMJgNmSvr/875xz0Tv8HBGRJcU5\nR2OklMZIKZ/d0UgqZbx7fjjYEfzR26f5wzd6yMtxbFlT5e0I1lazZU0V+YvcaO5Oj/hfBvab2e87\n5w4CfwY8DYwDXzGzzpusFwX+XEf8IpKtxqeS7O9N8EaXNz5wuG+QlEFJQS47YmF2tVTz9KaV1Ifu\nbv7AggzuOucKgHPARjO76Jw7AvwM+DLQhjcO0GSzvOFcEr9z7nngeYCGhobtvb29c4pLRGSpMTMO\nnhnk9356gr85fjlYvnl1Jf/rxd139Z4Ldanno3hH+xf9x33An/qJvsM5lwKqgcs3e4NbMbOXgJfA\nO+K/m/cQEclEyZRx7PwwHd1xOrrj7O2N0z/qNZKLlBbQFg3TFgvzlD+pbKHdSeJ/Dvj+jMc/Ap4A\nXnPOrQMKgP55jE1EZEkan0py6MwgnT1xOnoS7O9NMDrh3VGsPlTMo2traIuFaY+F01L6OafE75wr\nAT4CfGnG4peBl/1LPpPA583MnHN1wHfN7Bl/3e8DjwHVzrk+4HfM7L/N4+8gIpJWw+NT7OtN0Nkd\np7MnzqEzQ0wmU4BX0fPprXW0Rb1Ev6qyOM3RzjHxm9kYELlh2STwq7O89hzwzIzHz91jjCIiGeXy\nyIR3NO8n+mPnh0kZ5OU4Nq2u5Au7orRFw7Q2hjKyx4+atImI3IKZcSZ+lY6eOB3dA3T2JOjuvwJA\nUX4O2xpCvPjEWtpjYbY2VC2Jm8FnfoQiIosolTJ+fmkkGIjt7IlzcXgCgMrifNqiIZ5tW0NbLMym\nukoK8pbezd6V+EUkq01Opzh8dojOnjid3XH29iYYuurNuF1ZUUR7LEJ7NER7LMLa2rIlfbOWa5T4\nRSSrjE1Os793kA4/0R84k2B8yhuIbaou5aObVgYDsfWh4kWvuFkMSvwisqwlrkx6R/N+aeWRs96t\nF3McbKir4Ln2BtqjYVqjYWrKC9Md7qJQ4heRZeXc4NWg4qajO84J/9aKBXk5bKmv4rc+1ERbNMz2\nxhDlRflpjjY9lPhFZMny7ph1JRiE7eiOc3bwKgBlhXlsbwzx6a2raY+F2by6kqL8zLprVroo8YvI\nkjGdTPGu3/qgsyfO3p4EA1e81gfVZV7rg994JEZbNMz9qyrSflesTKXELyIZa3wqycEzg3R2x+no\nibO/N8GVySTg3Qz9sfW1tMdCtEXDxNLQ+mCpUuIXkYwxdHWK/b0Jf7JUnMN911sf3LeynF/cVu/1\nuImGWVmZWXe1WkqU+EUkbS4NjwdllR09Cd67MIz5rQ8211fy69daH0RDVJVkXuuDpUqJX0QWhZnR\nOzAWJPrOnjg9A2MAFOfnsq2xii8/6bc+WBOiuEADsQtFiV9EFkQqZbx3YcSvn/eS/aURr/VBVUk+\nrY1hPrujkbZYmI11FYt++8FspsQvIvPCa30wSEd3IpgwNTLu9aBfVVnEw82RYEZsS83yaH2wVCnx\ni8hduTIxzf7TXg/6t7vjHDwzyMS0NxDbXFPKxx9Y5d1ZKrp8Wx8sVUr8IjInA6MTdPZcP5o/em44\naH2wsa6Sz+5opD0WojUaprosO1ofLFVK/CIyq77EmD8b1kv2J2e2PlhTxW9/qJn2WJhtjSHKCpVK\nlhL93xIRzIyTl0ZnVNwkgtYH5YV5bI+G+MVtq2mPhtlcX0lhnipuljIlfpEsNJ1McfTc8PtuH5gY\n83rQV5cVsiMW5jcfidEWC3PfSrU+WG6U+EWywPhUkgOnB4Mkv/90gjG/9UFjpIQn719BezRMWyxM\nNFKigdhlTolfZBkaGptib+91G3UHAAANUElEQVT1+vnDZ4eYShrOwfoV5fyf2+uD0soVFWp9kG2U\n+EWWgYvD4+9rTXz84ghmkJ/r2Ly6ki/ujrEjFmZ7Q5jKkuzsQS/XKfGLLDFmRs/AWNCxsqM7zum4\n1/qgpCCXbQ0hntns1dBvWVOl1gfyAUr8IhkumTLeuzAcJPrOngSX/dYHoZJ82qJhPvdwI23RMBvU\n+kDmQIlfJMNMTCc53DfE2/6lm309CUYmvNYHq6uK2dUcCVoTN6v1gdwFJX6RNBudmGZfbyI4oj94\nZpBJv/VBS20ZH3+wLrjZSH2oJM3RynKgxC+yyPpHJ9g7Y0bs0XNDpAxycxwb6yr4tYcaaY+FaW0M\nEVHrA1kASvwiC8jM6EtcDfrbvN0d59TlKwAU+q0PXni8hbZYmK0Nan0gi0Nbmcg8SqWMk5dH6ei+\nPiP2/NA4AOVFebRFw3xm+xraYyE2rVbrA0mP2yZ+59x64AczFjUB3zCz7zjnXgReAKaBV8zsq7Os\n/zTwu0Au8F0z+9a8RC6SAab81gcd3QN0dCfY2xtn0G99UFteGAzCtkXDrF9ZrtYHkhFum/jN7Diw\nBcA5lwucBX7onHsc+BTwgJlNOOdqb1zXf/0fAB8B+oBO59yPzezdefwdRBbN1ckkB04n/LLKOPt7\nB7k65bU+iEZK+Mj9K4Jk36jWB5Kh7vRSz5NAl5n1Oue+DXzLzCYAzOzSLK9vB06a2SkA59wf4+0s\nlPhlSRgcm2Sv34P+7e44R84OMZ3yWh/ct7KCX2qtpz0WoS0aolatD2SJuNPE/yzwff/ndcAjzrlv\nAuPAV8ys84bXrwbOzHjcB+yY7Y2dc88DzwM0NDTcYVgi8+PC0Lg/G3aAzu4Exy+OAF7rgwfqq/jN\nR5toj3o96CuL1fpAlqY5J37nXAHwSeBrM9YNAQ8BbcCfOOeazMxmrjbLW9ksyzCzl4CXAFpbW2d9\njch8MjO6+694A7H+pZszca8HfWlBLtsaQ97tA2Ne64OifA3EyvJwJ0f8HwX2m9lF/3Ef8Kd+ou9w\nzqWAauDyjHX6gDUzHtcD5+4hXpG7lkwZx84PB9U2nT1x+kcnAQiXFtAWDfH5h6O0x8JsWFVBnlof\nyDJ1J4n/Oa5f5gH4EfAE8Jpzbh1QAPTfsE4nsNY5F8MbFH4W+JW7D1dk7sankrzTNxR0rNzf+/7W\nB4+srfFbE4dorinTQKxkjTklfudcCV5lzpdmLH4ZeNk5dwSYBD5vZuacq8Mr23zGzKadcy8Ar+KV\nc75sZkfn91cQ8YyMT3mtD/xEf6hvKGh9sLa2jE9sqWNHzCutrKsqTnO0Iunj3n9JPjO0trba3r17\n0x2GZLjLI17rg2vNzI6dHw5aH2yqq6DdT/Kt0TDh0oJ0hyuyoJxz+8ysdS6v1cxdWRKutT54u/va\nzcDjnOr3Wh8U5eewdU2IF55YS3s0zNaGKkrV+kDkpvTXIRkplTJ+fmnE71jpda68MOy1PqjwWx/8\nUtsa2qJhNq+upCBPA7Eic6XELxlhKpni8Nmh4Gi+syfB0FWv9cGKisLg/rDtsTDrasvVg17kHijx\nS1qMTU5z4PRg0MzswJkE41PeQGysupSnN64MWh+sCRer4kZkHinxy6JIXJkMauc7ehIcndH6YMOq\nCp5ta/B60EdD1Jar9YHIQlLilwVxbvBqUFbZ2RPn5xdHASjIzeHBNZU8/2gTbbEw2xtDVBSp9YHI\nYlLil3tmZnRdvuId0fvtD/oSXuuDssI8tjWG+OSDdbRFwzyo1gciaafEL3dsOpni2PkRr7+Nf0Q/\ncMVrfRApLaAtGubXd8XYEQtz38pytT4QyTBK/HJb41NJDp0ZDJqZ7e9NcGXS60FfHyrmQ+trvJuN\nxMI0VZdqIFYkwynxywcMj0+xrycRHNG/0zfEZNKruFm/opxf2LY6KK9cVanWByJLjRK/cGlknM7u\n6z1ujl0YxgzychybVlfyhV1Rr/VBY4iQWh+ILHlK/FnGzDgdH5vRmjhB94zWB9saQvzjJ9bSHvNa\nH5QUaBMRWW70V73MpVLG8YsjwdF8R3ecSyMTAFQW59MWDfFs2xraY2E2ra4kXwOxIsueEv8yMznt\ntT64dkS/tyfO8LjXg35lRRE7miJe64NomLW1ZWp9IJKFlPiXuCsT0+w/nQjq5w+eGQxaHzTVlPLM\n5lXBQGx9SK0PRESJf8mamE7yS//1LY6cHSKZMnIcbKir4Ln2Btr9HvQ15YXpDlNEMpAS/xJVmJdL\nU3Upu1sitEW91gflan0gInOgxL+E/adf3pLuEERkCVIJh4hIllHiFxHJMkr8IiJZRolfRCTLKPGL\niGQZJX4RkSyjxC8ikmWU+EVEsowzs3TH8AHOuctAb7rjmKEa6E93EDfIxJggM+PKxJggM+PKxJgg\nM+PKtJgazaxmLi/MyMSfaZxze82sNd1xzJSJMUFmxpWJMUFmxpWJMUFmxpWJMc2VLvWIiGQZJX4R\nkSyjxD83L6U7gFlkYkyQmXFlYkyQmXFlYkyQmXFlYkxzomv8IiJZRkf8IiJZJqsSv3PuM865o865\nlHOu9YbnHnDOvek/f9g5V+Qv/yvn3CF/+X9xzuXO8r7OOfd7zrmTzrl3nHPbZjz3eefcCf/f5+81\nJudciXPuFefce/7yb93kd/2sc+7gjH8p59wW/7nXnHPHZzxXu4hxRZ1zV2d89n+Z8dx2//1O+t+n\nu2HdhYrpI865ff56+5xzT8x4Lm3flb/+1/zv47hz7h/MWP60v+ykc+5f3GTdu9nev+mcO+OcG71F\nTHe9bS1gTHe9XS1wXPe0bS0YM8uaf8D9wHrgNaB1xvI84B3gQf9xBMj1f67w/+uA/x94dpb3fQb4\nS/81DwFv+8vDwCn/vyH/59C9xASUAI/7ywqAvwc+epvfezNwasbj933WfHxXc40LiAJHbvKZHcDD\n/vf4lzeuv4AxbQXq/J83AWcz5LvaABwCCoEY0OWvn+v/3OSvfwjYME/b+0PAKmB0jn9Td7RtLVRM\n97JdLXBc97RtLdS/rLoDl5kdA2a74fhTwDtmdsh/3cCMdYb9H/Pw/shmGxT5FPA98/5vvuWcq3LO\nrQIeA/7azOL+5/418DTw/XuIaQz4G3/ZpHNuP1B/m1/9uZmfOReLFFfA/74qzOxN//H3gE/j/aEu\naExmdmDGw6NAkXOu0Mwm5hL7An5XnwL+2I+j2zl3Emj3nztpZqf8z/1j/7Xv3mNcmNlbN1nnZu5o\n21qkmAJz2a4WMq573bYWSlZd6rmFdYA55151zu13zn115pPOuVeBS8AI8D9nWX81cGbG4z5/2c2W\n33NMflxVwCeAn97mvX6ZD/5x/qF/evmvZzv1XeC4Ys65A865v3XOPeIvW433/VyTru/q/wAO3PCH\nma7vaiG2qznFdQfma9uaj5jme7uar7iumc9t654suyN+59xPgJWzPPV1M/uzm6yWB+wG2vCOxn7q\nnNtnZj8FMLN/4F/X+/+AJ4C/vvFjZ3lPu7Z8Rkw1QMo594V7jck5l4f3B/d71478ZuOc2wGMmdmR\nGYs/a2ZnnXN/A3wV+JJzbnDG8wsZ13mgwcwGnHPbgR855zbywe+qBKh2zl2LezG+q43Af8A7yrsm\nnd/Vzbar4IDN/77WAyXOucfmI665utm2Bfx3oI4Pfl8LGdMttys/3p/gXR6buV0tdFzXPvtW21Y5\n3mXkXwO+dyfve7eWXeI3sw/fxWp9wN+aWT+Ac+4vgG3MOAozs3Hn3I/xTqlvTPx9wJoZj+uBc/7y\nx67F5Jz7r8BrZjaXU+PbxfQScMLMvnOb93mWG47IzOys/9/H/Z1Qq5m9MIeY7jku/2hnwv95n3Ou\nC++oqg+oN7P7/Pd9Du+7+9JCx+SvUw/8EPicmXXNiDdt3xU33664ttzMPuyc+5r/87+fp7jm6mbb\n1rXt/QvM/fu6p5hut135yz98h9vVPcflr3O7bWvEOfdHeJfxFiXx61KP51XgAedVW+QBHwLedc6V\n+dcIrx2dPQO8N8v6PwY+5zwPAUNmdt5/36eccyHnXAhvb//qvcTkx/LvgErgn9zqDZxzOcBngD+e\nsSzPOVft/5wPfBw4Mvs7zH9czrka51dGOeeagLV4g4PngRHn3EP+Ke/ngJsdhc13TFXAK8DXzOyN\nGcvT+l3hbVfPOucKnXMxvO+qA+gE1jrnYs65ArwE/OP5iGuuFmDbuqeYFmi7mo+4Fmrbujd3Ohq8\nlP8Bv4C3B58ALgKvznjuV/EGX44A/9FftgLvj+wd/7n/DOT5z/0W8Fv+zw74A7xKi8O8vyrgi8BJ\n/9+vz0NM9Xin+8eAg/6/3/Cf+yTwb2es/xjw1g2fVwrsm/E7/S5+lcJixIV3nfMoXiXKfuATM963\n1X/PLuD38ScYLkJM/wq4MuM1B4HadH9X/uOv+9/HcWZUo+AdhPzcf+7r87G9+8v/o79Oyv/vv5nP\nbWuhYrqX7WqB47qnbWuh/mnmrohIltGlHhGRLKPELyKSZZT4RUSyjBK/iEiWUeIXEckySvwiIllG\niV9EJMso8YuIZJn/Dd9rafcN7FTYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xf27b61c5f8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.plot(cur_data.variables['Longitude_corners'][:][0] , cur_data.variables['Latitude_corners'][:][1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#filter out all lon corners and lat corners in Illinois\n",
    "#then find OCO2 footprints\n",
    "\n",
    "\n",
    "# oco_data = Dataset('05\\\\oco2_LtSIF_160504_B7302r_160720173141s.nc4')\n",
    "#then find related fpars\n",
    "#then find related cdls\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "solar_zenith_angle\n",
      "sensor_zenith_angle\n",
      "sensor_azimuth_angle\n",
      "solar_azimuth_angle\n",
      "time\n",
      "longitude\n",
      "latitude\n",
      "footprint\n",
      "sounding_id\n",
      "measurement_mode\n",
      "surface_altitude\n",
      "daily_correction_factor\n",
      "reduced_chi2_757nm\n",
      "reduced_chi2_771nm\n",
      "continuum_radiance_757nm\n",
      "continuum_radiance_771nm\n",
      "SIF_757nm\n",
      "SIF_771nm\n",
      "uncorrected_SIF_757nm\n",
      "uncorrected_SIF_771nm\n",
      "SIF_757nm_uncert\n",
      "SIF_771nm_uncert\n",
      "SIF_757nm_relative\n",
      "SIF_771nm_relative\n",
      "uncorrected_SIF_757nm_relative\n",
      "uncorrected_SIF_771nm_relative\n",
      "orbit_number\n",
      "IGBP_index\n",
      "footprint_vertex_latitude\n",
      "footprint_vertex_longitude\n",
      "[6 7 8 ..., 3 4 6]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([ 12936.,  15732.,  15046.,      0.,  14882.,  15806.,      0.,\n",
       "         16179.,  16275.,  16585.]),\n",
       " array([ 1. ,  1.7,  2.4,  3.1,  3.8,  4.5,  5.2,  5.9,  6.6,  7.3,  8. ]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAD8CAYAAACcjGjIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAFcRJREFUeJzt3X2wXHd93/H3pxYG7IRI4As1khyJ\nRnFrmDQ4inHClEntYMvAIP8BM3IerKGa0Qw1FJJ2wG7+cAu4g9sMpp6AM6qtIKfUwjEk1gSBUY0p\nZQY/yDb4Sbi6sYl1sUGXkW0gNBDBt3/sT81GZ++Ddq+191rv18ydPed7fmfPd+9I+7nnYfekqpAk\nqd8/GncDkqTFx3CQJHUYDpKkDsNBktRhOEiSOgwHSVKH4SBJ6jAcJEkdhoMkqWPZuBsY1mmnnVZr\n1qwZdxuStKTce++9362qibnGLdlwWLNmDXv37h13G5K0pCT56/mM87CSJKnDcJAkdRgOkqQOw0GS\n1GE4SJI6DAdJUofhIEnqMBwkSR2GgySpY8l+QlqSxmnN5Z8dy3a/+eE3H5ftuOcgSeowHCRJHXOG\nQ5LtSQ4meeio+ruTPJrk4ST/ua9+RZLJtuzCvvqGVptMcnlffW2Su5LsT/KpJCcv1IuTJA1nPucc\nPgH8EXDjkUKSfwlsBH6pqn6U5OWtfhawCXg18Ergfyb5xbbax4A3AlPAPUl2VdUjwNXANVW1M8kf\nA1uA6xbixUl6fhvXcf8TwZzhUFVfTrLmqPI7gQ9X1Y/amIOtvhHY2eqPJ5kEzmnLJqvqMYAkO4GN\nSfYB5wG/1cbsAP4DhoO0pPgm/fwz7DmHXwT+RTsc9L+S/GqrrwQO9I2barWZ6i8Dnqmqw0fVJUlj\nNOylrMuAFcC5wK8CNyd5FZABY4vBIVSzjB8oyVZgK8AZZ5xxjC1LkuZr2D2HKeAz1XM38FPgtFZf\n3TduFfDkLPXvAsuTLDuqPlBVbauq9VW1fmJizrvcSZKGNOyew1/QO1fwpXbC+WR6b/S7gP+R5CP0\nTkivA+6mt4ewLsla4Fv0Tlr/VlVVkjuAtwE7gc3ArSO8Hmnsxnn8/Xh9QErPf3OGQ5KbgN8ATksy\nBVwJbAe2t8tbfwxsrqoCHk5yM/AIcBi4rKp+0p7nXcBtwEnA9qp6uG3i/cDOJB8C7gduWMDXt6j4\npiFpqZjP1UqXzLDod2YYfxVw1YD6bmD3gPpj/P0VTZKkRcBPSEuSOgwHSVKH4SBJ6jAcJEkdhoMk\nqcNwkCR1GA6SpA5vE3qCeL7f0lDSwnLPQZLU4Z6Dnre8x4A0PPccJEkdhoMkqcNwkCR1GA6SpA7D\nQZLUYThIkjrmDIck25McbHd9O3rZv0tSSU5r80lybZLJJA8kObtv7OYk+9vP5r76ryR5sK1zbZIs\n1IuTJA1nPnsOnwA2HF1Mshp4I/BEX/kieveNXgdsBa5rY19K7/air6N317crk6xo61zXxh5Zr7Mt\nSdLxNWc4VNWXgUMDFl0DvA+ovtpG4MbquRNYnuR04EJgT1UdqqqngT3AhrbsJVX11XYP6huBi0d7\nSZKkUQ11ziHJW4FvVdXXj1q0EjjQNz/VarPVpwbUJUljdMxfn5HkFOAPgAsGLR5QqyHqM217K71D\nUJxxxhlz9ipJGs4wew7/BFgLfD3JN4FVwH1J/jG9v/xX941dBTw5R33VgPpAVbWtqtZX1fqJiYkh\nWpckzccxh0NVPVhVL6+qNVW1ht4b/NlV9W1gF3Bpu2rpXODZqnoKuA24IMmKdiL6AuC2tuz7Sc5t\nVyldCty6QK9NkjSk+VzKehPwVeDMJFNJtswyfDfwGDAJ/DfgXwNU1SHgg8A97ecDrQbwTuD6ts5f\nAZ8b7qVIkhbKnOccquqSOZav6Zsu4LIZxm0Htg+o7wVeM1cfkqTjx09IS5I6DAdJUofhIEnqMBwk\nSR2GgySpw3CQJHUYDpKkDsNBktRxzF+893yw5vLPjrsFSVrU3HOQJHUYDpKkDsNBktRhOEiSOgwH\nSVKH4SBJ6jAcJEkd87kT3PYkB5M81Ff7L0m+keSBJH+eZHnfsiuSTCZ5NMmFffUNrTaZ5PK++tok\ndyXZn+RTSU5eyBcoSTp289lz+ASw4ajaHuA1VfVLwP8BrgBIchawCXh1W+fjSU5KchLwMeAi4Czg\nkjYW4GrgmqpaBzwNzHYbUknScTBnOFTVl4FDR9W+UFWH2+ydwKo2vRHYWVU/qqrH6d0X+pz2M1lV\nj1XVj4GdwMYkAc4Dbmnr7wAuHvE1SZJGtBDnHP4V8Lk2vRI40LdsqtVmqr8MeKYvaI7UJUljNFI4\nJPkD4DDwySOlAcNqiPpM29uaZG+SvdPT08fariRpnoYOhySbgbcAv11VR97Qp4DVfcNWAU/OUv8u\nsDzJsqPqA1XVtqpaX1XrJyYmhm1dkjSHocIhyQbg/cBbq+qHfYt2AZuSvDDJWmAdcDdwD7CuXZl0\nMr2T1rtaqNwBvK2tvxm4dbiXIklaKPO5lPUm4KvAmUmmkmwB/gj4WWBPkq8l+WOAqnoYuBl4BPg8\ncFlV/aSdU3gXcBuwD7i5jYVeyPx+kkl65yBuWNBXKEk6ZnPez6GqLhlQnvENvKquAq4aUN8N7B5Q\nf4ze1UySpEXCT0hLkjoMB0lSh+EgSeowHCRJHYaDJKnDcJAkdRgOkqQOw0GS1GE4SJI6DAdJUofh\nIEnqMBwkSR2GgySpw3CQJHUYDpKkDsNBktQxnzvBbU9yMMlDfbWXJtmTZH97XNHqSXJtkskkDyQ5\nu2+dzW38/nb/6SP1X0nyYFvn2iRZ6BcpSTo289lz+ASw4aja5cDtVbUOuL3NA1xE777R64CtwHXQ\nCxPgSuB19O76duWRQGljtvatd/S2JEnH2ZzhUFVfBg4dVd4I7GjTO4CL++o3Vs+dwPIkpwMXAnuq\n6lBVPQ3sATa0ZS+pqq9WVQE39j2XJGlMhj3n8IqqegqgPb681VcCB/rGTbXabPWpAXVJ0hgt9Anp\nQecLaoj64CdPtibZm2Tv9PT0kC1KkuYybDh8px0Soj0ebPUpYHXfuFXAk3PUVw2oD1RV26pqfVWt\nn5iYGLJ1SdJchg2HXcCRK442A7f21S9tVy2dCzzbDjvdBlyQZEU7EX0BcFtb9v0k57arlC7tey5J\n0pgsm2tAkpuA3wBOSzJF76qjDwM3J9kCPAG8vQ3fDbwJmAR+CLwDoKoOJfkgcE8b94GqOnKS+530\nroh6MfC59iNJGqM5w6GqLplh0fkDxhZw2QzPsx3YPqC+F3jNXH1Iko4fPyEtSeowHCRJHYaDJKnD\ncJAkdRgOkqQOw0GS1GE4SJI6DAdJUofhIEnqMBwkSR2GgySpw3CQJHUYDpKkDsNBktRhOEiSOgwH\nSVLHSOGQ5PeSPJzkoSQ3JXlRkrVJ7kqyP8mnkpzcxr6wzU+25Wv6nueKVn80yYWjvSRJ0qiGDock\nK4F/A6yvqtcAJwGbgKuBa6pqHfA0sKWtsgV4uqp+AbimjSPJWW29VwMbgI8nOWnYviRJoxv1sNIy\n4MVJlgGnAE8B5wG3tOU7gIvb9MY2T1t+fpK0+s6q+lFVPU7v/tPnjNiXJGkEQ4dDVX0L+EPgCXqh\n8CxwL/BMVR1uw6aAlW16JXCgrXu4jX9Zf33AOpKkMRjlsNIKen/1rwVeCZwKXDRgaB1ZZYZlM9UH\nbXNrkr1J9k5PTx9705KkeRnlsNJvAo9X1XRV/R3wGeDXgeXtMBPAKuDJNj0FrAZoy38OONRfH7DO\nP1BV26pqfVWtn5iYGKF1SdJsRgmHJ4Bzk5zSzh2cDzwC3AG8rY3ZDNzapne1edryL1ZVtfqmdjXT\nWmAdcPcIfUmSRrRs7iGDVdVdSW4B7gMOA/cD24DPAjuTfKjVbmir3AD8aZJJensMm9rzPJzkZnrB\nchi4rKp+MmxfkqTRDR0OAFV1JXDlUeXHGHC1UVX9LfD2GZ7nKuCqUXqRJC0cPyEtSeowHCRJHYaD\nJKnDcJAkdRgOkqQOw0GS1GE4SJI6DAdJUofhIEnqMBwkSR2GgySpw3CQJHUYDpKkDsNBktRhOEiS\nOgwHSVLHSOGQZHmSW5J8I8m+JL+W5KVJ9iTZ3x5XtLFJcm2SySQPJDm773k2t/H7k2yeeYuSpONh\n1D2H/wp8vqr+KfDPgX3A5cDtVbUOuL3NA1xE7/7Q64CtwHUASV5K725yr6N3B7krjwSKJGk8hg6H\nJC8B3kC7R3RV/biqngE2AjvasB3AxW16I3Bj9dwJLE9yOnAhsKeqDlXV08AeYMOwfUmSRjfKnsOr\ngGngT5Lcn+T6JKcCr6iqpwDa48vb+JXAgb71p1ptprokaUxGCYdlwNnAdVX1WuBv+PtDSINkQK1m\nqXefINmaZG+SvdPT08faryRpnkYJhylgqqruavO30AuL77TDRbTHg33jV/etvwp4cpZ6R1Vtq6r1\nVbV+YmJihNYlSbMZOhyq6tvAgSRnttL5wCPALuDIFUebgVvb9C7g0nbV0rnAs+2w023ABUlWtBPR\nF7SaJGlMlo24/ruBTyY5GXgMeAe9wLk5yRbgCeDtbexu4E3AJPDDNpaqOpTkg8A9bdwHqurQiH1J\nkkYwUjhU1deA9QMWnT9gbAGXzfA824Hto/QiSVo4fkJaktRhOEiSOgwHSVKH4SBJ6jAcJEkdhoMk\nqcNwkCR1GA6SpA7DQZLUYThIkjoMB0lSh+EgSeowHCRJHYaDJKnDcJAkdRgOkqSOkcMhyUlJ7k/y\nl21+bZK7kuxP8ql2lziSvLDNT7bla/qe44pWfzTJhaP2JEkazULsObwH2Nc3fzVwTVWtA54GtrT6\nFuDpqvoF4Jo2jiRnAZuAVwMbgI8nOWkB+pIkDWmkcEiyCngzcH2bD3AecEsbsgO4uE1vbPO05ee3\n8RuBnVX1o6p6nN49ps8ZpS9J0mhG3XP4KPA+4Kdt/mXAM1V1uM1PASvb9ErgAEBb/mwb///rA9aR\nJI3B0OGQ5C3Awaq6t788YGjNsWy2dY7e5tYke5PsnZ6ePqZ+JUnzN8qew+uBtyb5JrCT3uGkjwLL\nkyxrY1YBT7bpKWA1QFv+c8Ch/vqAdf6BqtpWVeurav3ExMQIrUuSZjN0OFTVFVW1qqrW0Duh/MWq\n+m3gDuBtbdhm4NY2vavN05Z/saqq1Te1q5nWAuuAu4ftS5I0umVzDzlm7wd2JvkQcD9wQ6vfAPxp\nkkl6ewybAKrq4SQ3A48Ah4HLquonz0FfkqR5WpBwqKovAV9q048x4Gqjqvpb4O0zrH8VcNVC9CJJ\nGp2fkJYkdRgOkqQOw0GS1GE4SJI6DAdJUofhIEnqMBwkSR2GgySpw3CQJHUYDpKkDsNBktRhOEiS\nOgwHSVKH4SBJ6jAcJEkdhoMkqWPocEiyOskdSfYleTjJe1r9pUn2JNnfHle0epJcm2QyyQNJzu57\nrs1t/P4km2fapiTp+Bhlz+Ew8G+r6p8B5wKXJTkLuBy4varWAbe3eYCL6N0feh2wFbgOemECXAm8\njt4d5K48EiiSpPEYOhyq6qmquq9Nfx/YB6wENgI72rAdwMVteiNwY/XcCSxPcjpwIbCnqg5V1dPA\nHmDDsH1Jkka3IOcckqwBXgvcBbyiqp6CXoAAL2/DVgIH+labarWZ6pKkMRk5HJL8DPBp4L1V9b3Z\nhg6o1Sz1QdvammRvkr3T09PH3qwkaV5GCockL6AXDJ+sqs+08nfa4SLa48FWnwJW962+CnhylnpH\nVW2rqvVVtX5iYmKU1iVJsxjlaqUANwD7quojfYt2AUeuONoM3NpXv7RdtXQu8Gw77HQbcEGSFe1E\n9AWtJkkak2UjrPt64HeBB5N8rdX+PfBh4OYkW4AngLe3ZbuBNwGTwA+BdwBU1aEkHwTuaeM+UFWH\nRuhLkjSiocOhqr7C4PMFAOcPGF/AZTM813Zg+7C9SJIWlp+QliR1GA6SpA7DQZLUYThIkjoMB0lS\nh+EgSeowHCRJHYaDJKnDcJAkdRgOkqQOw0GS1GE4SJI6DAdJUofhIEnqMBwkSR2GgySpY9GEQ5IN\nSR5NMpnk8nH3I0knskURDklOAj4GXAScBVyS5KzxdiVJJ65FEQ7AOcBkVT1WVT8GdgIbx9yTJJ2w\nFks4rAQO9M1PtZokaQyWjbuBJgNq1RmUbAW2ttkfJHl0yO2dBnx3yHWPt6XUKxzVb64eYydzW9K/\n20EW0e/7efe7XSxy9ci9/vx8Bi2WcJgCVvfNrwKePHpQVW0Dto26sSR7q2r9qM9zPCylXmFp9buU\neoWl1e9S6hWWVr/Hq9fFcljpHmBdkrVJTgY2AbvG3JMknbAWxZ5DVR1O8i7gNuAkYHtVPTzmtiTp\nhLUowgGgqnYDu4/T5kY+NHUcLaVeYWn1u5R6haXV71LqFZZWv8el11R1zvtKkk5wi+WcgyRpETmh\nwiHJ9iQHkzw07l7mkmR1kjuS7EvycJL3jLunmSR5UZK7k3y99fofx93TfCQ5Kcn9Sf5y3L3MJsk3\nkzyY5GtJ9o67n7kkWZ7kliTfaP9+f23cPQ2S5Mz2Oz3y870k7x13X7NJ8nvt/9hDSW5K8qLnbFsn\n0mGlJG8AfgDcWFWvGXc/s0lyOnB6Vd2X5GeBe4GLq+qRMbfWkSTAqVX1gyQvAL4CvKeq7hxza7NK\n8vvAeuAlVfWWcfczkyTfBNZX1dK4Dj/ZAfzvqrq+XX14SlU9M+6+ZtO+wudbwOuq6q/H3c8gSVbS\n+791VlX93yQ3A7ur6hPPxfZOqD2HqvoycGjcfcxHVT1VVfe16e8D+1iknxqvnh+02Re0n0X9V0eS\nVcCbgevH3cvzSZKXAG8AbgCoqh8v9mBozgf+arEGQ59lwIuTLANOYcDnwRbKCRUOS1WSNcBrgbvG\n28nM2iGarwEHgT1VtWh7bT4KvA/46bgbmYcCvpDk3vYtAYvZq4Bp4E/aIbvrk5w67qbmYRNw07ib\nmE1VfQv4Q+AJ4Cng2ar6wnO1PcNhkUvyM8CngfdW1ffG3c9MquonVfXL9D7dfk6SRXvYLslbgINV\nde+4e5mn11fV2fS+tfiydnh0sVoGnA1cV1WvBf4GWNRfwd8Ofb0V+LNx9zKbJCvofSHpWuCVwKlJ\nfue52p7hsIi14/efBj5ZVZ8Zdz/z0Q4hfAnYMOZWZvN64K3tWP5O4Lwk/328Lc2sqp5sjweBP6f3\nLcaL1RQw1bfneAu9sFjMLgLuq6rvjLuROfwm8HhVTVfV3wGfAX79udqY4bBItZO8NwD7quoj4+5n\nNkkmkixv0y+m94/4G+PtamZVdUVVraqqNfQOJ3yxqp6zv8BGkeTUdkEC7fDMBcCivdquqr4NHEhy\nZiudDyy6iyiOcgmL/JBS8wRwbpJT2vvD+fTORT4nTqhwSHIT8FXgzCRTSbaMu6dZvB74XXp/1R65\n1O5N425qBqcDdyR5gN73ZO2pqkV9eegS8grgK0m+DtwNfLaqPj/mnubybuCT7d/DLwP/acz9zCjJ\nKcAb6f0Vvqi1vbFbgPuAB+m9fz9nn5Y+oS5llSTNzwm15yBJmh/DQZLUYThIkjoMB0lSh+EgSeow\nHCRJHYaDJKnDcJAkdfw/4tVvNohBBf0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xf27b8f41d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i in oco_data.variables:\n",
    "    print(i)\n",
    "    \n",
    "    \n",
    "# print(oco_data.variables['footprint_vertex_latitude'][:].shape)\n",
    "# print(oco_data.variables['footprint_vertex_longitude'][:].shape)\n",
    "print(oco_data.variables['footprint'][:])\n",
    "\n",
    "plt.hist(oco_data.variables['footprint'][:])\n",
    "\n",
    "# plt.plot(oco_data.variables['footprint_vertex_longitude'][:][0], oco_data.variables['footprint_vertex_latitude'][:][0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def find_inner_regions(areas1, areas2):\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'FparExtra_QC': (('YDim:MOD_Grid_MOD15A2H', 'XDim:MOD_Grid_MOD15A2H'),\n",
       "  (3678, 10145),\n",
       "  21,\n",
       "  3),\n",
       " 'FparLai_QC': (('YDim:MOD_Grid_MOD15A2H', 'XDim:MOD_Grid_MOD15A2H'),\n",
       "  (3678, 10145),\n",
       "  21,\n",
       "  2),\n",
       " 'FparStdDev_500m': (('YDim:MOD_Grid_MOD15A2H', 'XDim:MOD_Grid_MOD15A2H'),\n",
       "  (3678, 10145),\n",
       "  21,\n",
       "  4),\n",
       " 'Fpar_500m': (('YDim:MOD_Grid_MOD15A2H', 'XDim:MOD_Grid_MOD15A2H'),\n",
       "  (3678, 10145),\n",
       "  21,\n",
       "  0),\n",
       " 'LaiStdDev_500m': (('YDim:MOD_Grid_MOD15A2H', 'XDim:MOD_Grid_MOD15A2H'),\n",
       "  (3678, 10145),\n",
       "  21,\n",
       "  5),\n",
       " 'Lai_500m': (('YDim:MOD_Grid_MOD15A2H', 'XDim:MOD_Grid_MOD15A2H'),\n",
       "  (3678, 10145),\n",
       "  21,\n",
       "  1)}"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "raw_data.datasets()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from shapely.geometry import Point\n",
    "from shapely.geometry import Polygon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "mat= '48.6666045 −0.2645971422 −0.2546445294 −0.1379386595 \\\n",
    "−0.2645971 0.0022624549 0.0010580484 −0.0001827183 \\\n",
    "−0.2546445 0.0010580484 0.0014918324 0.0008762569 \\\n",
    "−0.1379387 −0.0001827183 0.0008762569 0.0037290135'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'48.6666045 −0.2645971422 −0.2546445294 −0.1379386595 −0.2645971 0.0022624549 0.0010580484 −0.0001827183 −0.2546445 0.0010580484 0.0014918324 0.0008762569 −0.1379387 −0.0001827183 0.0008762569 0.0037290135'"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "mat = mat.replace('−', '-')\n",
    "m = mat.split()\n",
    "m = np.array(m)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "m = m.astype(float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "m = m.reshape((4,4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "e = '5.943 \\\n",
    "−0.143 \\\n",
    "−7.040 \\\n",
    "−3.729 \\\n",
    "3.282 \\\n",
    "1.130 \\\n",
    "0.557'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "e = e.replace('−', '-')\n",
    "e = e.split()\n",
    "e = np.array(e)\n",
    "e = e.astype(float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "y = '10 \\\n",
    "5 \\\n",
    "−6 \\\n",
    "25 \\\n",
    "20 \\\n",
    "31 \\\n",
    "5'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "y = y.replace('−', '-')\n",
    "y = y.split()\n",
    "y = np.array(y)\n",
    "y = y.astype(float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "beta = '10.72 \\\n",
    "0.15 \\\n",
    "0.10 \\\n",
    "−1.56'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "beta = beta.replace('−', '-')\n",
    "beta = beta.split()\n",
    "beta = np.array(beta)\n",
    "beta = beta.astype(float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "111.16541199999999"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SSE = np.sum(e**2)\n",
    "SSE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>\n",
    "a)</p>\n",
    "\n",
    "<table>\n",
    "<tr>\n",
    "<td>Source</td>\n",
    "<td>Sum of Sqaures</td>\n",
    "<td>df</td>\n",
    "<td>Mean</td>\n",
    "<td>Squares F</td>\n",
    "</tr>\n",
    "\n",
    "<tr>\n",
    "<td>Regression</td>\n",
    "<td>903.69173</td>\n",
    "<td>3</td>\n",
    "<td>301.2306</td>\n",
    "<td>8.129</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>Error</td>\n",
    "<td>111.16541</td>\n",
    "<td>3</td>\n",
    "<td>37.0551</td>\n",
    "<td></td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>Total</td>\n",
    "<td>1014.8571</td>\n",
    "<td>6</td>\n",
    "<td></td>\n",
    "<td></td>\n",
    "</tr>\n",
    "\n",
    "\n",
    "</table>\n",
    "\n",
    "\n",
    "\n",
    "<p>b) <br>\n",
    "By doing the previous F test in a), we have obatined a test statistics of 8.129. Since the critcal value at p = 0.05 is 9.2766. We do not reject the null hypothesis, and there is no evidence that at least one estimator has a non-zero coefficient</p> \n",
    "\n",
    "\n",
    "<p>c)<br>\n",
    "The unbiased variance for y is 37.05. The variance for estimated $\\beta_1$ is 0.08383557703. Hence the test statistics is 0.518056443932. The degree of freedom is 3. The critical t value is 3.182. Hence we do not reject the null hypothesis, and there is no evidence that $\\beta_1$ is not zero\n",
    "</p>\n",
    "\n",
    "\n",
    "<p>d)<br>\n",
    "The variance for estimated $\\beta_3$ is 0.1381791073. The test statistics for $\\beta_3$ is 4.19665634218 which is greater than the critical value of 3.182. Hence we reject the null hypothesis and conclude that $\\beta_3$ is not 0. \n",
    "\n",
    "</p>\n",
    "\n",
    "\n",
    "<p>e)<br>\n",
    "The model explains the variation quite well since the $R^2$ is 0.8904. \n",
    "\n",
    "</p>\n",
    "\n",
    "\n",
    "<p>f)<br>\n",
    "\n",
    "We are going to use F test to compare the reduced model and full model. \n",
    "$$F_{test} = \\frac{MS_{diff}}{MSE_{full}}$$\n",
    "\n",
    "We have $SS_{diff} = 121.6974 - 111.16541 = 10.5319\\\\\n",
    "         MS_{diff} = \\frac{10.5319}{(4-2)} = 5.266$\n",
    "         \n",
    "<br>\n",
    "Hence F statistics is 5.266/37.0551 = 0.14211 with degree of freedom (4-2, 3) which is (2,3). And we do not reject the null hypothesis. Hence we conclude that there is no evidence that the full model explains the variation better than the reduced model. Hence that person's argument is correct. </p>\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.14211269164028703"
      ]
     },
     "execution_count": 104,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SST = np.sum((y - np.mean(y))**2)\n",
    "(121.6974 - 111.16541)/2\n",
    "5.266/37.0551"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "SSR = np.sum((y - np.mean(y))**2) - SSE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df_r = len(y) - len(beta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df_t = len(y) - 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "301.230576952381"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MSR = SSR/df_r\n",
    "MSR\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "37.055137333333327"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MSE = SSE/df_r\n",
    "MSE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.1292527468628748"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MSR/MSE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "37.0551373333\n"
     ]
    }
   ],
   "source": [
    "sigma_2 = np.sum((y - np.mean(y))**2) / (len(y) - 4)\n",
    "\n",
    "\n",
    "sigma_2 = np.sum(e**2)/(len(e) -4)\n",
    "print(sigma_2)\n",
    "var_beta_0 = m[0,0]*sigma_2\n",
    "var_beta_1 = m[1,1]*sigma_2\n",
    "var_beta_2 = m[2,2]*sigma_2\n",
    "var_beta_3 = m[3,3]*sigma_2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1803.34771329\n",
      "0.08383557703\n",
      "0.0552800544603\n",
      "0.13817910736\n"
     ]
    }
   ],
   "source": [
    "print(var_beta_0)\n",
    "print(var_beta_1)\n",
    "print(var_beta_2)\n",
    "print(var_beta_3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.19665634218\n"
     ]
    }
   ],
   "source": [
    "sd_beta_3 = np.sqrt(var_beta_3)\n",
    "test_stats = 1.56 /sd_beta_3\n",
    "\n",
    "print(test_stats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 10.,   5.,  -6.,  25.,  20.,  31.,   5.])"
      ]
     },
     "execution_count": 91,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  4.86666045e+01,  -2.64597142e-01,  -2.54644529e-01,\n",
       "         -1.37938660e-01],\n",
       "       [ -2.64597100e-01,   2.26245490e-03,   1.05804840e-03,\n",
       "         -1.82718300e-04],\n",
       "       [ -2.54644500e-01,   1.05804840e-03,   1.49183240e-03,\n",
       "          8.76256900e-04],\n",
       "       [ -1.37938700e-01,  -1.82718300e-04,   8.76256900e-04,\n",
       "          3.72901350e-03]])"
      ]
     },
     "execution_count": 92,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.89046200957207211"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SSR/(SSR+SSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}