Match-up_tool_marine_sciences.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: right;\" src=\"./img/es2_logo_tiny.jpg\"></img>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1 align=\"center\">Match Up Tool - Marine Sciences</h1>\n",
    "<h2 align=\"center\">Example Generating regression plot of insitu data vs remote sensed</h2>\n",
    "<hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"./index.ipynb\"><< Index </a>&nbsp;<space>&nbsp;<space>&nbsp;<space>&nbsp;<space>&nbsp;<space>&nbsp;<space>&nbsp;<space>&nbsp;<space>&nbsp;<space>&nbsp;<space> <a href=\"./Point-extract_example_marine_sciences.ipynb\">02-Single-point extraction - Marine Sciences >></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Background"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>Match up services are used by scientists to both validate remote sensed data as well as in-situ data. Match ups are generated from a file of latitude, longitude and date points. These can often hold several hundred points depending on the length and purpose of the cruise that generated the in-situ data. Current methods for finding EO data matchups normally take the form of custom written scripts that have to find the most appropriate file, and within it the most appropriate pixel. This can be quite lengthy operation even when you discount the time it takes to create the script.</p>\n",
    "<p>By providing a tool that can lookup values from the Ocean Colour CCI datasets we can provide a service to reduce the time needed to produce matched up datasets</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Addressed Community\n",
    "\n",
    " * Marine and Earth Science researcher\n",
    " \n",
    "### Benefits\n",
    "\n",
    "<p>Web Coverage Service 2.0 with its processing extension Web Coverage Processing Service (WCPS) offers an efficient retrieval of single points of data without the need to download extraneous data.</p>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Produce a CSV file including the matchups from the OC-CCI daily dataset\n",
    "\n",
    "### Description\n",
    "\n",
    "<p> A researcher has a Comma Seperated Values (CSV) file that contains locations, dates and insitu chlorophyll measurements. They would like to query the OC-CCI Daily data to generate a corosponding CSV file that has the EO data points added</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>Data required</h3>\n",
    "<div class=\"alert alert-info\">\n",
    "<ul type=\"square\">\n",
    "<li><a href=\"\">OC CCI Daily </a> -  Chlorophyll values matching the dates in CSV</li>\n",
    "<li><b>Format:</b> NETCDF</li>\n",
    "<li><b>Spatial resolution:</b> Lat/Lon grid, 4Km x 4Km</li>\n",
    "<li><b>Temporal resolution:</b> Daily</li>\n",
    "</ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3> Workflow - traditional </h3>\n",
    "<div class=\"alert alert-success\">\n",
    "<ul type=\"square\">\n",
    "<li>(1) Retrieve data that corresponds to your spatial and temporal area</li>\n",
    "<li>(2) Using a scripting language find which file is needed for  each data point</li>\n",
    "<li>(3) extract pixel value using scripting language</li>\n",
    "<li>(5) Plot the match up regression</li>\n",
    "</ul>\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This has the major overhead of needing to download a very large dataset just to get access to individual pixels. Depending on the spatial and temporal spread of inistu data points the data download could be several hundred Gb or even Tb's.\n",
    "\n",
    "Using WCPS combined with python to parse the results can remove this need and speed up the match up process dramatically. The mechanism to do this is outline below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "#### Load required libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "from datetime import datetime, timedelta\n",
    "\n",
    "import csv\n",
    "import requests"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define our  WCPS query\n",
    "\n",
    "<p>We will use a single WCPS query multiple times, the query template will hold the coverage we want to extract from placeholders for Latitude, Lonitude and Date</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "query_template = '''\n",
    "for c in (V3_daily_cci_all_insitu_match_test)\n",
    "return\n",
    "encode (\n",
    "   c[Lat({}), Long({}), ansi(\"{}\")], \"csv\"\n",
    ")'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#### Define some helper functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# This will simple create a new WCPS query using our template and given values\n",
    "def create_query(lat,lon,date):\n",
    "   \"\"\"simply tests POST wcps requets\"\"\"\n",
    "   query = query_template.format(lat,lon,date)\n",
    "   return query\n",
    "\n",
    "# This will send the WCPS query off to the web service and clean the result to a single value\n",
    "def run_query(lat,lon,date):\n",
    "   \"\"\"simply tests POST wcps requets\"\"\"\n",
    "   query = query_template.format(lat,lon,date)\n",
    "   payload = {'query' : query}\n",
    "   r = requests.post('http://earthserver.pml.ac.uk/rasdaman/ows/wcps'\\\n",
    "   , data=payload)\n",
    "   return r.text[1:-1]\n",
    "# This is just used to clean the dates from the CSV file and make them the format required by the web service\n",
    "def clean_date(date):\n",
    "   \"\"\"simply return string in rasdaman approved format\"\"\"\n",
    "   return datetime.strptime(date, \"%d/%m/%Y\").strftime(\"%Y-%m-%d\")\n",
    "\n",
    "# Main bulk of work is here, we loop through CSV file and get value then add back to new CSV file\n",
    "# if the value is null we do not add it\n",
    "def loop_csv(csv_file):\n",
    "   \"\"\"loops through csv list and tests each line\"\"\"\n",
    "   with open(csv_file, \"r\") as csvfile:\n",
    "      with open('output_'+csv_file, 'w') as csv_out:\n",
    "         csvFile = csv.reader(csvfile, delimiter=',')\n",
    "         csv_writer = csv.writer(csv_out)\n",
    "         head = next(csvFile)\n",
    "         head.append(\"match_value\")\n",
    "         csv_writer.writerow(head)\n",
    "         x = 1\n",
    "         for row in csvFile:\n",
    "            d = clean_date(row[2].split(' ')[0])\n",
    "            #d = clean_date(row[2])\n",
    "            x = x + 1\n",
    "            #q = {'query' : create_query(row[0],row[1],d)}\n",
    "            #queries.append(q)\n",
    "            match_value = run_query(row[0],row[1],d)\n",
    "            \n",
    "            if(match_value != '9.96921e+36'):\n",
    "                #print(match_value)\n",
    "                row.append(match_value)\n",
    "                csv_writer.writerow(row)\n",
    "                \n",
    "                \n",
    "                \n",
    "                \n",
    "def trendline(xd, yd, order=1, c='r', alpha=1, Rval=False):\n",
    "    \"\"\"Make a line of best fit\"\"\"\n",
    "\n",
    "    #Calculate trendline\n",
    "    coeffs = np.polyfit(xd, yd, order)\n",
    "\n",
    "    intercept = coeffs[-1]\n",
    "    slope = coeffs[-2]\n",
    "    power = coeffs[0] if order == 2 else 0\n",
    "\n",
    "    minxd = np.min(xd)\n",
    "    maxxd = np.max(xd)\n",
    "\n",
    "    xl = np.array([minxd, maxxd])\n",
    "    yl = power * xl ** 2 + slope * xl + intercept\n",
    "\n",
    "    #Plot trendline\n",
    "    #plt.plot(xl, yl, c, alpha=alpha)\n",
    "\n",
    "    #Calculate R Squared\n",
    "    p = np.poly1d(coeffs)\n",
    "\n",
    "    ybar = np.sum(yd) / len(yd)\n",
    "    ssreg = np.sum((p(xd) - ybar) ** 2)\n",
    "    sstot = np.sum((yd - ybar) ** 2)\n",
    "    Rsqr = ssreg / sstot\n",
    "\n",
    "    if not Rval:\n",
    "        #Plot R^2 value\n",
    "        plt.text(0.8 * maxxd + 0.2 * minxd, 0.8 * np.max(yd) + 0.2 * np.min(yd),\n",
    "                 '$R^2 = %0.2f$' % Rsqr)\n",
    "    else:\n",
    "        #Return the R^2 value:\n",
    "        return Rsqr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To run our match up now we just need to call the `loop csv` method with a path to a CSV file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "loop_csv('cleaned_cci_db_4.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Producing the match up regression plot\n",
    "<p>To do this we will use the data we generated in the previous step with a library called matplotlib.</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.871923697539\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x4a37690>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh0AAAGHCAYAAAAUdeBfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3Xd4VFX6wPHvJCSQRkIJHRNBumJIAClKE4EgCLiABlC6\nBaX+RNdGsSysCouirBRBFIyAsKg0AQWkCgRC7x1pARJKgADJ+f1xZibTkkyGSSYJ7+d55kluf++d\nW957zrl3QAghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBC\nCCFEwTMKSPN0EFkoBHwCnAJSgf95NhyXrQZWuThtGjDSoruXsd8D9xZSjvkWOObpIPKgUeT9482d\nVuP6Pi+c1wy9XzXxcBwu8/J0ALmgF/pLMn3uAH8D35F3T+Q5QRk/eVkf4A1gLvAiMN6JaToBS4EE\nIAX93c4BmjsYtzTwGbAfSAauA1uBd4Fgi/FWA7tcWQGje93Wef17spXf4s0NOX28dQMGZ2P848Cv\nORMKYL++5dCJ16M5sKw0YGIGwzqTdy7KO4ETWYyzHjjH/XEtvm/0Qu+E76IP1D7AVHTycRzw91Rg\nucwb8PV0EFn4ETjp5LgGYAb6u90K/BP9Xb8DbDH2b2gxfj10YpIMTAZeMn6mAteA3yzGXY0+Ybhq\nNfCHi9P6or8rk15ISUd+lNPH2yKyt92PAb/kUCygSykLWXTXRe+3L+bAstKALzIYlpeSjjfRsTyR\nwfBwdInuhGzMsxl5Z/1cUijrUQqMpcA24//TgYvAW0BH4IdcjMMAFAZu5eIyQe/cqbm8zOwqBVxx\nctz/A3oC/zH+b+lfQHfgrrE7BF1VcweoAxy0GHcKOiHt51rIbnfb0wEIt8iN4y07JSmGHItCu5tB\n/5xebl7gB9x00P8HYAz6Znetg+Ex6O0zO+dCy3vu5yKddca/FW36Vwd+Ai6hd6QtQHsH09cG1gA3\n0G0Q3gV6Y39XehxdrNkafUd+A3jZOCwEneWeQichh9DZse2B+jwQB1xFX5R3AoMshvug2wEcMsZ8\nEb2Tt7QYZxT2dcyFgPeBI8blHwM+xv4OzbQOjwObjcs4AryAcwKAcRbruR/rRCHcGFszoBbpVWEZ\nZfN+wNvAPnR1jCOz0d8d6O1dDhiGdcJhcgGdqLjiJfS2uAH8heO7Gh/gA/R3mISu1vkTvb62bNt0\n2JqJLrFxdMOwHL1tM/IlulTHz8GwWOAs6fteB2AxurrqFnAYeI+szxnNcPzdhRv797Tp78zx5sz+\n7UgxdHXaLvR6XwGWoI9dW2HokoBk4Dy6aq+1g3V5ApiHLja/hS6ZGw8UsZnfKOyPN1O1QEdgt3H6\n3cblWApCnxeOG8c5j/5u6xiHrwbakr5N08h+aZNp2v8jfR++hT6+69qMWwZdqnjaOM4ZYCF6m5ms\nJr1NRzPjfCC9NNKy1OO4sb8ty3m4UxVgPnr/vok+D8UCRW3G64E+Rm+g98dYoIKDGHcBUehjOJmM\nzx2njeN0xvHx2g19XG1Bb8tJwAHj8i+iq5nDHExn6zjOb8/CwGjjck3777+xP+c/hb5GJqKPnf3o\na8M9u59KOmyFG/+es+hXC13HdgqdoSYDz6EPsH8Y/wKUR3+Zqegd7gb6Tvk29ncfCqiGznq/Rhft\nH0BX66wByhr7nQQaG5dbFhhqnP4p47Qr0VUBADWBRqQXMY5CVy9MRR/swegTRx3jdJaxWJqGPhHM\nAz4FGqAv5jWAZ22me8g43jT0Dt4XXbQeB+wlYwb0ybyZcdp4oI1xeeXRicAFdALzLjpBeds4bUYX\n0MfRF5TxDtbJkWfQ39FPToybHX3R3+l6YyyVgZ+By1hXEwUbx/0B/V0XNXb/BtQHdtjMN7N1+g69\nrVqjkwKTMuh2LKMymfZHYADwNNbbwh99oZ9useye6CR3HDpJehKdOBVFJ8auslw3Z4+3UTi3f9uq\njE6e5qIvymXQCega9DF01jheALo6rDT6Yn8efUFo4WCendEJxiT0hekxYCD64tQ1k3U1eRx9bH2F\n3q6D0BfEB9D7Deh96h/oBGUvUBJ9bqgObAc+Qje4rgAMMU5zPZPtkJlu6CTnv8buN4EFQCXSSy/m\no7fXF+gLXGl0wleR9DYLlm069gIj0PvLZNLv8jc4GNdSTrSD8UUfZz7o+M+ht9vT6P3oqnG8d43x\nzkGXfpZCf69/ovczUwmsAkqgk9dY9PF4PpPlzzbOz/Z4fQS9/482dtdDVwf/gE5WHgReRScONXFc\nkmLi7Pb0Qp+LG6O/l33oBHwoUBXdPg5jXIvQ5+r30W3lqqCvOcIJvdBZdgv0wVsBfUBfQGfsQRbj\nrkRvaB+beaxDJwomX6ATDss7pmLok1Aq9iUdaejkwdJ76Ayysk3/f6GrAcobuyegs83Miinjybq+\ndhTWd16PGrsn24z3CemlDibHjf0aW/QriT4QPs1iuR2M075t038ueltVsui3GufaUgwyzvMZJ8YF\nfTLfluVY2YvDB32yicM6ee9njM2yTYcX9gl+MPqiN82mfxr6hG3SC+vSMwM6oYm1mW4oentmdWd0\nCp08WuqC/fdre+cO+sJ0Hevj41us77KbkXlJh2Udv7PHmzP7tyOO2lSEoffb9yz6DTPGZlnCUhh9\n8bRdF0fb5S30trcsNR2F45KOm+gLiskjxv6vWfRLIuM2CyaLgKNZjGPpONbbMNy43AtYN6Jub+z/\ntLE7xNg9LIv5r8Z6n8+sTccxdIKb1Twykp02HRHG7mczGB/0PnEX/T1aqoW+kbQ8d602zq+/E3GC\n3n43sa9CGWOcz0PGbkf71WPGcXpY9GuG/T7p7PbsgV5P2+ThJazbwA0xdhd3MM97dj9Vr6xEH2An\n0SfdU+ii0mvG4cXRd4rz0AdhSYvPcnSmV9Y4bht01m55YUoEZuE4OTgKrLDp1wWdRSfZLOt3dCM0\n006VCAQCrTJZt0TgYdJ3YGe0Nf61fUJknPHv0zb996DvSk0uoi8MD5K5tugd3fYkMQ69raKdCdaG\nqVj0WqZjWY/v7LjOqguEou9KLeuzv8W+XUqaxThe6H3NB52wRGZzuQp9AnsGvV+YdEd/P1m1lp+H\n/k4CLPo9h767svx+LdscBaH3zXXoUpHq2YzZkewcb67s32DdPsYbfYeajK5is9zubdDrb/l0Rwrp\nJYuWLLdLgDHejeh9OcKJmFZinaTtQt9tWx5HiehSx7LkvDlY76+mamdTPDfR27E5+gKa35jWrQ2O\nqxVBJyQGdOmf5X54Hl0NYfsk3C0cV2c4koQuFXmG9IcWDOgq8y3G+ZvmaeKD3lePGKevg3t0QZdu\nHMB6PS2rxUDvf6CrAd2eI9xPSccAdJFgZ/ROEIEu2jZ5CL0zfIhOTiw/o9An+1LGccNI31ksHclg\n2Y7qW6ugL7gJNstaYbOsSeiT5FJ0ovQN9nXAI9AnhIPoROgT9B1UZsLQF0Pb9TiP3tFtn5Zw9FRJ\nErqEJ6vlnEGf7C2Zqk5ceSrDVCQalOlY1uM7O66zTCUKh2z638XxHWhP9HdjapNwAX3xt61XdsZ3\n6BOoqTi0Gvoi+r0T084xTmsqJQo0xmFb+lEL3fg2CX3ivmAx/2DuXXaON1f2b4zzH4r+jm6Rfqw9\ngvV2D8Pxseuo3wPoxPIyOpG9gL6jBOe2i6PjKBHr4+hNdJJ1Ct1OaCRZJ/euso3HdMExxZOCLgGI\nRp8b1gDD0VUseZmpWuEY+saqH/q4W4a+Flh+/1XQ+8oh7PfF6uibC0t/k3HDWUdmoxPUDsbuRuh9\nzrL0ww9dvWNq92baV0Nwz/EGej1rYX/NOYD18TYHfQMyDV0dFYtOWNzSKPh+atOxmfQi9oXojH4S\nur7vMukJ2KdYPz5pyXQSym69o6P6OAP6ju6TDKYxXcwS0AlSa/SBH41usPoduugddJ2pqf66FfoA\nGwq8gk5SMuNoXRztXBm1xHdmR3R3C/Z9xr+1ca7YfT+6OskHXXWV02zXtwf6zuh/6EZbF0ivcqpE\n9u1Dl5L0QCcCPdAXh7lOTPsXuqi9K/pk0h5dtDvHYpwQ9MUlCeuGxlHG+DO7Wcno2PC26c7O8ebq\n/m2qp/8GncxfNsY3AdduuLyN8wlBF4+b3vdSAZ2IODNPZ46jeeh17oRe3+HoC/+z6IumOzkTz+fo\nUqCO6PPQh+h9twW66iu7MttHnHmhWgoZl1qYShMsSw7eQH8/pv3nC3T8DdAJhJcxpjY43h627WUy\na1/hyCJ04t4Nfcx1QyctP1qMMxF9Pv8PuuTMVELzI1nvV85uTy900p5RVdkp499b6JL25ugS7zbo\n0tA/0NvvfnrpnUt6oTeSbTF2U2P/D4zdpYzdzrTQPUB6MaSlL3D89IqjC+NurIuznWVA162nkfEF\nKwB9UTpl0W8U1jvL28Zu26Ly0sb+lsnQcRyvw2qyroOdjL7QB9r0N9VXDrCZnzNtOvzQ7Wf24NyJ\n/p/GZT3vxLjOxtHQOM+XbPr7oO8WLbfLQuxLREB//7alIlm16TAZiN6uZdAX5+w0kh2LPnEGGWOz\nvaPvaFzm4zb9+2Nfn/wt1iV5tXHc3qYF1nX82TnebDnavx2Jx3FD09NYfz/LM5iXqa2HbfuAHjbj\nPYV9+4VROG7T4agtQkZ18iahxvgsH7v8Ffe06XB0AbLdB209hL4QW5asrcZ6m0aRcZuOOBy/bfgk\nzrXp2I/jR1BB709ppN+1O2I6dj80dg83dldxYtmrce0dPlPRyVIp9E2HbfKYhH37riLo5MRy32iG\n/THo7PZcRNbHTEZM14snXZze7H6qXrG1Bl36MQD95ZqKSV9Gn8htWRax/YbecS3ftlccXa/ubCnI\nXOM8HLXVCCH9ztC2MY8i/W2ZhY1/S9iMk4y+kGT2ciJTS+ohNv2H2Qy/V4vR6/K6Tf+h6J14qQvz\nvIm+465h/OtID3SLcNDtLs6i25E4OrGUwrphoTO2oEuhXsG6IWQv7ItDU9HJouXd42NYv7wsu2LR\n+8IX6KL3WdmYdg563+mJvnO1LSEx3e1Znh98sU4QLVnu8yeM0ze1Gcd22uwcb67s36BP2LbnuC7o\nx6ctLUM33LZMlIpg31jQ0XYxkL03g2bFC/v9JwG9/1qubzK508bCD/tGjkfRSUdm299Uneqo+vUI\nupTB8rhph/3jqRlZYpze9kYyBH0O3o7ev0An1rYl+rvR5x5T/PPR362jR9UNuKdB5Wz0+k5Bt6Ow\nbVjqaF8d6KCfI85uz7no/dxRI1g/0kuJHH1npifs7vmFd/dT9Yojn6KLMvuiH2F7DV2CsQudmR5D\n3/k3RH9ZpoZin6AvaivQxWI3jPM4gf7CnEk8PkWf5Bah7xa3oe/gHkE/XROGLg7+xjjPP9BFgWHo\nnXE76dUMe9GNgbYZp6lL+iN3GdmJfufDS+iD9U90G5cX0VnzGifWAbKuOvnVGNvH6LurnehE6xl0\nUaJtexdnq2I+RddP/h+6GPAndJ1zGfSduukRNNB3EZ3QJ6t49AXaVNUWiS4BMT3O52wcd9GJymT0\ndzMXffHvhT4pW07/q3H5/zPG8CD6YrsH+xIgZ5nqpzujS1aykyRuR7fl+Rc6+ZhjM3y9cZ4zSb8z\nz+ydLJbregV9TA1EHwdH0SdA23pxcP54c2X/Bn1sjUDfKW5EH1vdsP9+JqOT4lh0VcI59MXLVIxu\nOp73oU/wnxnju2aM414v/paxFEWXxMxDHyvX0W3R6mJdKrEVXUU2zvj/dXLmNefV0I3b56DX/y56\nXw7FunoArNfD1AjyFWNsycAmdInLNPR+uwy9npXR2/sIzh3/Y0lviG96BUE59LFXGut3wTyJfj/N\nXHRpYyH0vnwXnWyA3h/eQ1eZhaMfe7+GPk47ohOFcRbzdKW6eA36ezU9vr/AZvgiY1xX0Nu5oTH2\nS04sz9nt+T16n/kafc7cgL4hrI7enq3Qx9hI9EMWi9GlJaXQNw2ncFzCL2z0Qmexjp4SMDUeOkz6\nl/MgOgk4gy4OO4neCTvZTPsoeke6aRznHfSJNg3rE2xmrx8OQF+MD6Lr0S6giw2Hkp4Qmupxz5H+\nAq9JWBcfvoM+oC+jD+496CoFy3r0kdjXV3qTXmefgj4hfIT9I4wZrcMqnCsONb0c7LRxOftxXKy7\niuwXXZq2z0V0K/u/0ScY2ztt0AnJOOPyb6BPhnHoun/LhqbZieMV9Pa7iW4v0RjH2+Wf6O14E32R\niEa383CmesX2MWwT0+OB/3UwLCsfGqc9kMHwhuiTUjLp79F4yhiLZdGuo3UogT75XUd/L5PQ7xpw\nVNzuzPHmzP7tiC86Of3bOJ0psXb0/YSjL9qml4ONQ+9baaSXmIE+QS9HN06+gD6Bmx57tVw3R8eb\nM9UrPujSu+3oC9A19IXgZZtp/NHJ82XjfLOqarE9hsNxrnqlOOnvC7mGTkY3oJMtS462aXt0qcJt\n9Law3D5D0fvVTdLfheHs+QR0kjHFOI/b6NKgn7H+rkCv5zT0ed700q2VOP5tpk7GWK4ZP3vQ35fl\nU1OunKNM/o3eDraPu4Mu3foGvU9dRd+cVMW+6q0Z9scgOL89C6Grk3YZx72ELvF/j/RzYHP0DZLp\nZXCn0Pua7esdRB4wAX3Suh9e/Ss8z/QOlMZZjShcYnpfQW48uiqEyCVN0HcYf6MP8A6Zjw7oIuGP\nSX9N8DH0Ex25ybb1dAl01phRS3wh3G0RjhuoiuyzPZ6LoIu5M3utvBAim/JCmw5/dFHiN+h6Lmfa\nQ8xFV2H0QVeNlCXrolZ324guvtqPrkfsi66f/zCziYRwg+fRT4m0xfo3eITrFqDbZO1AF3X3QBdv\nd/dkUEKInOXoUTtbbdD1ip5+O97H6PrwZHTd9Roc/1aDEO6Whq7vn8L9/QSaOw0m/YfhbqCfTuri\n0YiEEDnOmaRjEvqpkbHohi4H0I3FHL27XgghhBB5RF6oXsmuSuiXFt1EP84Uik5ESqCrW4QQQgiR\nB+XHpMMLXSLSnfQf8RqGfk/Dq+jH7myVRVqgCyGEEK44a/zcs/yYdJxFP9Nv+auh+9GPqVbA/pXO\nZcuVK3fmzJkzuRSeEEIIUaD8jX4Hyj0nHvkx6ViHfilSAOmv2q2KLv047WD8smfOnGHWrFnUqFEj\nl0LMWUOGDGHChAmeDsMtCtK6gKxPXlaQ1gVkffKygrQu+/bto0ePHuXRtQUFIukIwPr3MCqhX398\nifQ3IZYj/dW2P6DfojkD/da/UHRD0m9wXLUCQI0aNYiMdPRS0vwnJCRE1iWPkvXJuwrSuoCsT15W\nkNbF3fLC43b10K/53YZ+R8d44/+jjcPLABUtxk9Gv445BP066Vno19/K+wqEEEIUOHdS7/DJ+k9I\nupXk6VDuWV4o6VhN5smPozeNHsDxr7MKIYQQBcr55PN8uuFTqpWoRofqzry0O+/KC0mHEMKBQ4cO\nce3ataxHzEFJSUls27Yt6xHzgYK0LiDrk5e5Y12UUiw+tJjm4c0J8A1g/hPzCbwRmGPbKCgoiCpV\nqmQ94j2SpCMfiomJ8XQIblOQ1gXctz6HDh2iatWqbpnXvYqKivJ0CG5TkNYFZH3yMnety0hGumU+\nzjh48GCOJx6SdORDBelCXZDWBdy3PqYSjoL01JUQIm8yPqGSKyWrknQIkYcVpKeuhBAiLzy9IoQQ\nQoj7gCQdQgghhMgVknQIIYQQIldI0iGEKJBGjRpFnTp1PB2Gy8LDw/n88889HYYQbiVJhxAiz+vV\nqxedOnXydBi5ymAwYDAYPB2GEG4lSYcQQgghcoUkHUIIt2nWrBmDBg1iyJAhFC9enDJlyjB9+nSu\nX79Oz549KVq0KFWqVGHZsmXmadLS0ujbty+VKlXC39+f6tWr88UXX5iHjxo1iu+++46ff/4ZLy8v\nvLy8+PPPPwE4ffo0MTExlChRgsDAQOrVq8fmzZutYvr+++8JDw8nJCSEmJgYrl+/bh7mqAojIiKC\n0aNHm7u9vLz4+uuviY6Oxt/fn8qVKzN//vwMt8GUKVMoX748Simr/h06dKBv374AHDlyhA4dOlCm\nTBmCgoKoX78+v//+e4bzPH78OF5eXuzcudPcLykpyWpbAOzevZvo6GiCgoIoU6YML774IpcuXTIP\n/+mnn3jkkUfw9/enZMmSPPXUU9y4cSPD5QrhbpJ0CJFP9ekD1arpv3lpfjNnzqRUqVJs2bKFgQMH\n8vLLL9O5c2eaNGnC9u3badWqFS+88AI3b94EdNJRsWJFfvrpJ/bt28eIESN45513mDdvHgDDhw+n\na9euREdHc+7cOc6dO0fDhg25fv06TZs25ezZs/z666/s2rWLt99+m7S0NHMsR44c4ZdffmHJkiUs\nWrSINWvWMHbsWPNwR1UYjvq9//77dOnShZ07d9K9e3eef/559u/f73D9u3btyqVLl1i1apW53+XL\nl/ntt9/o0aMHAMnJybRr144//viD+Ph42rRpQ/v27Tl16pTL2z0pKYkWLVoQFRVFXFwcy5Yt4/z5\n83Tt2hWAs2fPEhMTQ79+/di/fz+rV6/mH//4h11yJIS4N5GAiouLU0LkF3FxcSqz/bZ3b6VCQ5UC\n/bd373tbnrvm17RpU9WkSRNzd2pqqgoMDFQ9e/Y09zt37pwyGAzqr7/+ynA+r732murcubO5u2fP\nnqpjx45W40yePFkVLVpUJSYmOpzHyJEjVUBAgLp+/bq535tvvqkaNGhg7g4PD1eff/651XQRERFq\n9OjR5m6DwaAGDBhgNU6DBg3s+lnq2LGj6tu3r1WsFSpUyHB8pZR6+OGH1ZdffukwtmPHjimDwaB2\n7NhhHp6YmKgMBoNas2aNUkqpDz/8ULVu3dpqnqdOnVIGg0EdOnRIxcXFKYPBoE6cOJFpHOL+k9n5\nxjTMeC29Z1LSIUQ+tH49JCTo/xMSdHdemJ/BYKB27drmbi8vL0qUKMEjjzxi7leqVCkALly4YO73\n1VdfERUVRalSpQgKCmLq1KlZ3vXHx8cTGRlJSEhIhuOEh4cTEBBg7i5TpozVcp3VsGFDu+59+/Zl\nOH737t2ZP38+d+7cAWD27Nk8//zz5uHXr1/njTfeoGbNmhQrVoygoCD27dt3TyUdO3bsYNWqVQQF\nBZk/NWrUwGAwcOTIESIiInjyySd55JFH6Nq1K9OmTSMpKf//VLrIXyTpECIfatwYQkP1/6Ghujuv\nzM/Hx8eq22AwWPUzVV2YqkF+/PFHhg8fTv/+/VmxYgU7duygd+/epKSk2M3Hkr+/f5ZVA45isax+\n8fLyspuHKVHITFbLbdeuHUopFi1axKlTp1i3bh3du3c3D3/jjTdYuHAhY8aMYd26dcTHx/PII49w\n+/Zth/Pz8vKyW65tnMnJyTzzzDPs2LHD6nPo0CGeeOIJvLy8WLFiBUuXLqVmzZpMnDiRatWqcfz4\n8SzXVwh3kaRDiHxo+nRo1w6qVtV/p0/PW/PLjvXr19OoUSNeeeUVHn30USpVqsThw4etkgxfX1/u\n3r1rNV3t2rWJj48nMTHR5WWHhoZy5swZc/fVq1c5duyY3XgbN2606t60aRM1a9bMcL5FihTh2Wef\nZfbs2cTGxlK9enUiIiLMwzds2EDv3r3p0KEDtWrVonTp0g6XaxknYBVrfHy81TiRkZHs3r2bsLAw\nKlWqZPXx9/c3j9eoUSNGjRrF9u3b8fX1ZeHChRkuVwh3k6RDiHxq+nQ4cMB9CYI75qeUynbDxKpV\nq7J161aWL1/OwYMHef/999m6davVfB588EF27tzJwYMHuXjxInfv3iUmJoYyZcrQsWNHNmzYwNGj\nR5k/fz6bNm1yetktWrTg+++/Z926dezatYuePXvi7e1tN95PP/3EjBkzOHjwICNHjmTr1q28/vrr\nmc67e/fuLFq0iBkzZliVcgBUqVKF+fPnm0sjunXrlul28/Pzo0GDBowdO5b9+/ezZs0a3nvvPatx\nXnvtNS5fvkxMTAxbt27lyJEj/Pbbb/Tp04e0tDT++usv/vWvfxEXF8fJkyeZP38+CQkJ8ivGIldJ\n0iGEcBtXXmj18ssv8+yzz/Lcc8/RoEEDEhMTGTBggNV8+vfvT7Vq1ahbty6lS5dmw4YN+Pj4sHz5\nckqVKkXbtm2pXbs2n3zyCYUKFcowFtt+b7/9Nk2bNqVdu3a0b9+eTp06UblyZbsYR48ezY8//sij\njz7KrFmzzKUXmWnRogUlSpTg4MGDdOvWzWrY+PHjKVasGI0aNaJDhw5ER0dn+WvC06dP5+7du0RF\nRTFs2DA+/vhjq3UpW7Ys69evJzU1lVatWlG7dm2GDh1KsWLF8PLyIjg4mLVr19K2bVuqVavGiBEj\nGD9+PK1bt850uUK40/3wurtIIC4uLk5+IlzkG9u2bTM/+ij7rWd5eXmxcOFCnnnmGU+HIkSOyOx8\nYxoGRAHb7nVZUtIhhBBCiFwhSYcQQgghckUhTwcghBB5meUjtkKIeyMlHUIIIYTIFZJ0CCGEECJX\nSNIhhBBCiFwhSYcQQgghcoUkHUIIIYTIFZJ0CCGEECJX5IWkownwK/A3kAZ0yMa0jYG7wPYciEsI\ncR/r1asXnTp18nQYQhQoeSHp8EcnDa8Zu539tagQ4DtgZTamEUIIp0ycOJGZM2fm+HK8vLzMn6Cg\nICIiIuyWu2DBAp566ilKlSpFcHAwjRo1Yvny5Tke28mTJ3n66acJCAigdOnSvPnmm6SmpmY6zZkz\nZ+jWrRtlypQhMDCQqKgoFixYYB6+evVqq3W2/MTFxZnHGzRoEHXr1qVw4cLUqVPH4bLmzp1LREQE\nAQEBhIeH89lnn1kNd2a7LViwgLp161KsWDECAwOpU6cOs2bNslvW33//TY8ePShZsiT+/v7Url3b\nKt7z58/Tq1cvypcvT0BAANHR0Rw+fNg8/Pjx4xmu9/z5883jhYeH2w3/5JNPzMMvX75MmzZtKF++\nPEWKFOGBBx5g4MCBXLt2LVvbRmhpgLM/cPAjMBoYSeYlHZGAiouLU0LkF3Fxcaog7LcpKSk5Nu/b\nt2/n2LyIi/80AAAgAElEQVRzk8FgUDNnzlTnz59XBw8eVEOGDFEGg0H9+eef5nGGDBmiPv30U7V1\n61Z1+PBh9c477yhfX1+1ffv2HIvr7t276uGHH1atWrVSO3bsUEuXLlWhoaHqnXfeyXS65s2bq8ce\ne0xt2bJFHTt2TH300UfK29vbHOvt27fV+fPnrT79+/dXlStXtprPoEGD1KRJk9SLL76o6tSpY7ec\nJUuWKB8fHzV58mR17NgxtXjxYlWuXDn15ZdfmsdxZrutXr1aLVy4UO3fv18dPXpUff7556pQoUJq\n2bJl5nEuX76swsLCVJ8+fdSWLVvU8ePH1YoVK9SRI0eUUkqlpaWpBg0aqKZNm6qtW7eqAwcOqJdf\nflmFhYWp5ORkpZRSqampduv9wQcfqKJFi5rHUUqp8PBw9dFHH1mNZzk8MTFR/fe//1VxcXHq5MmT\n6vfff1fVq1dXzz//fLa2jaXMzjemYcZraYHjbNLRG9iELqkZhSQdooDJr0lH06ZN1euvv64GDx6s\nSpYsqVq0aKGUUmrXrl2qTZs2KjAwUJUuXVq98MIL6uLFi+bprl69qrp166YCAgJU+fLl1RdffKGa\nNm2qhgwZYh4nLCxMffjhh+qFF15QRYsWVb1791ZKKbV27Vr1+OOPKz8/P1WxYkU1aNAgq5P0V199\npR566CFVpEgRVbp0adW5c2fzsHnz5qmHH35Y+fn5qRIlSqiWLVuap+3Zs6fq2LGjedxbt26pgQMH\nqlKlSqkiRYqoxx9/XG3ZssU8fNWqVcpgMKjff/9dRUVFKX9/f9WoUSN14MCBTLeZwWBQP//8s7n7\n2rVrymAwqHHjxmU6Xa1atdQHH3yQ6Tj3YsmSJcrb21tduHDB3O/rr79WwcHB6s6dOxlOFxgYqGbN\nmmXVr0SJEuqbb75xOP7t27dVaGio+uijjxwOHzlypIqIiLDrHxMTo7p27WrVb+LEiapixYoZxqaU\nc9stMjJSjRgxwtz91ltvqSZNmmQ4/oEDB5TBYFB79+4190tLS1OlSpVS06ZNy3C6iIgI1a9fP6t+\n4eHhasKECZnGZ+vzzz+3Wu/sbpvcTDryQvVKdlUBxgA90EmKECIPmTlzJkWKFGHDhg18/fXXJCUl\n0aJFC/OvWC5btozz58/TtWtX8zTDhg1j48aN/Prrr/z222+sXr2a+Ph4u5+m/+yzz6hTpw7x8fG8\n//77HDlyhOjoaLp06cKuXbuYM2cO69at4/XXXwdg69atDB48mI8++oiDBw+ybNkymjZtCsDZs2eJ\niYmhX79+7N+/n9WrV/OPf/wDpXRtrcFgsFr+m2++yYIFC/juu+/Ytm0bDz30EK1btyYxMdEqxvfe\ne4///Oc/bN26lUKFCtGnT58st5lpmXfv3mXatGkYDIZMf104LS2Na9euUaJEiUznGxgYSFBQUIaf\nAQMGZDjtxo0bqV27NqGhoeZ+rVq14urVq+zZsyfD6aKjo/nxxx9JTEwkLS2NH3/8kZSUFJo1a+Zw\n/F9++YXLly/Tu3fvTNfF1u3btylcuLBVvyJFinD69GlOnjzpcJqstptSit9//51Dhw7RokULqxij\noqLo0qULpUuXJjIykmnTppmHp6SkAFjFYzAY8PX1Zf369Q6XFRcXx44dO+jbt6/dsLFjx1KyZEki\nIyP57LPPMq3SOnPmDAsWLKBVq1bmfq5sm/tVViUd3sAW4GWLfqNwoqTjiSeeUO3bt7f6/PDDD9nK\nJoXILc6WdJy5ekbtPLfTrv/2s9vVuWvnrPolJCeouDP289tzYY86deXUvQVs1LRpUxUVFWXV78MP\nP1StW7e26nfq1CllMBjUoUOH1NWrV5Wvr6+aP3++efiVK1dUQECAGjp0qLlfWFiYevbZZ63m07dv\nX/Xyyy9b9Vu7dq3y9vZWKSkpav78+So4OFhdu3bNLta4uDhlMBjUiRMnHK6LZUnH9evXla+vr4qN\njTUPv3Pnjipfvrz69NNPlVLpJR1//PGHeZwlS5Yog8GQaTWTwWBQfn5+KjAwUHl7eyt/f381d+7c\nDMdXSql///vfqkSJEiohISHT8Y4cOZLpJ7Pp+/fvr9q0aWPVLzk5WRkMBquqB1s3btxQLVq0UAaD\nQfn4+Kjg4GC1YsWKDMePjo5WTz/9dIbDMyrpmDJligoICFC///67Sk1NVQcOHFDVq1dXBoNBbdq0\nyeG8MtpuSUlJKiAgQPn4+ChfX1+70onChQurIkWKqHfffVfFx8erKVOmKD8/PzVz5kyllC6tCQsL\nU127dlWJiYkqJSVFjR07VhkMBrttaPLqq6+qWrVq2fUfP368WrNmjdq1a5f6+uuvVbFixdSwYcPs\nxnv++eeVv7+/eRmW+1h2t43pfPPxxx/bXSefeOKJ+7p6JcQ4zh2LT6pFv2YOppHqFZHvOJt0jFw1\nUpUfV96uf9C/gtS4DdbF81PjpipGYTduza9qqqHLhtr1d0WzZs3USy+9ZNWvc+fOytfXVwUGBlp9\nvLy81LJly1R8fLwyGAzq1CnrxCcyMtIq6QgPD1f/+te/rMapW7euKly4sNV8AwIClJeXl9q/f7+6\ndu2aql27tgoNDVUvvPCCmj17trpx44ZSStext2zZUhUtWlR16dJFTZ06VSUmJprnbZl07NixQxkM\nBnXy5Emr5Xfq1En16dNHKZWedFhWG23bts3hulkyGAxq8uTJ6siRI2r58uXqoYceMs/TkdmzZ5sv\nKDmpf//+dsmiM0lHp06dVIMGDdQff/yhdu7cqUaPHq1CQkLUrl277MY9deqU8vb2VgsWLMhwfhkl\nHUrpag8/Pz9VqFAhVaJECTV69GhlMBjU5s2b7cbNbLulpaWpI0eOqB07dqhx48Ypf39/tXjxYvNw\nHx8f1bhxY6tpBg0apBo2bGjujouLUxEREcpgMKhChQqp6Oho1bZtW9W2bVu75d24cUMFBwer8ePH\nZ7jeJjNmzFA+Pj52bZjOnTunDhw4oH755RdVq1YtFRMT4/K2kTYdGTMANW0+XwH7jP/7O5hGkg6R\n7+TXko5mzZpZJQpK6TvZzp07O7zTTk5OzjDpqFOnjl3S8fnnn1uNU6NGDTV48GCH8zadpO/evatW\nrlyp3nzzTVW5cmVVpUoVlZSUZJ7H+vXr1ciRI1Xt2rVVqVKl1LFjx5RSziUdHTt2VH379lVKpScd\nV65cMQ/fvn17pqUpStm36di8ebMyGAxq/fr1duPGxsYqf39/tWTJkgznZykgIMAu2bP8vPrqqxlO\nO2LECLuL/dGjR5XBYFDx8fEOp9m7d68yGAxqz549Vv1btmypXnnlFbvxP/jgA1W6dGl19+7dDOPI\nLOlQSicMZ86cUXfu3DGXLFkmfkplf7v169dPtWzZ0twdFham+vfvbzXOpEmTVPny9gn/1atXzcuv\nX7++ev311+3G+e6775Svr69dnI7s3r1bGQwGdfDgwQzHWbdunTIYDOrvv/+26u/MtlEqd5OOvPDT\n9gHodhomlYAI4BJwCt1+oxzQE73ie22mTwBuOegvRIFXNqgsZYPK2vWPKBNh16+kf0lK+pe0618z\ntGaOxGYSGRnJ/PnzCQsLw9vb2254pUqV8PHxYfPmzVSoUAGAK1eucPDgwQzbAVjOe8+ePVSqVCnD\ncby9vXnyySd58sknGTlyJCEhIaxatYqOHTsC0KhRIxo1asSIESMICwtj4cKFDBkyxGoelStXxtfX\nl3Xr1hETEwPAnTt32LJlC8OGDcvO5shSvXr1aNWqFWPGjOHXX38194+NjaVv377MmTOH6Ohop+a1\nc+fOTIcXLVo0w2ENGzbk448/JiEhwdyuY8WKFQQHB1OzpuN9Ji1NN7Oz/Z69vLzM7VZMlFLMmDGD\nF1980eF+4SyDwUDZsvoYiI2NpVGjRlZtNlzZbqmpqeZ1AWjcuDH79++3GufgwYOEh4fbTRsUFATA\noUOHiIuL4+OPP7Yb55tvvqFDhw5ZtskBiI+Px8vLi1KlSmUar+Vfk6y2jSfkhaSjHvCH8X8FjDf+\n/y3QBygDVMxkelMWJoTwMKWU3cXltddeY+rUqcTExPDmm29SrFgxDh8+zJw5c/jmm28ICgqiZ8+e\nDB8+nOLFixMaGsrIkSPx9va2a0hq66233qJBgwYMHDiQvn37EhAQwN69e1m5ciUTJ05k0aJFHD16\nlCZNmlCsWDGWLFmCUopq1aqxefNmVq5cSevWrQkNDeWvv/4iISGBGjVq2C0nICCAV1991RxjxYoV\n+eSTT7h165bDhoD3asiQITz99NPs27ePGjVq8MMPP9CzZ0+++OIL6tWrx7lz5wDw9/fPNHHILBnL\nSqtWrahZsyYvvPACn3zyCWfPnuX999/ntddew8fHB9DvrnjyySf5/vvvqVevHtWrV6d69eq89NJL\nfPbZZxQvXpyFCxeycuVKFi9ebDX/P/74g+PHj9OvXz+Hyz98+DDXr1/n3Llz3Lx5kx07dqCUolat\nWvj4+HDp0iXmzZtHs2bNuHXrFjNmzGD+/PmsWbPGPA9nttuYMWOoV68elSpVIiUlhSVLljBr1iym\nTJlins/QoUNp1KgRY8aMoUuXLmzevJmpU6cydepU8zjz5s0jNDSUBx54gF27djF48GA6depEy5Yt\n7dZr7dq1LF261G6dN23axKZNm2jevDlBQUFs3LiRYcOG8cILLxAcHAzA0qVLOXfuHPXq1SMwMJA9\ne/YwfPhwmjRpQsWK+lLpzLYROUeqV0S+k18fmXVUvaKUUocOHVLPPvusKlasmPL391c1atSwahx3\n7do11b17dxUQEKDKlSunJkyYoB577DGrd0I4ql5RSqktW7aoVq1aqaCgIBUYGKgeffRRNWbMGKWU\nLnZu1qyZKl68uPL391cRERFq3rx5Siml9u3bp9q0aWN+BLZ69erqq6++Ms+3V69eqlOnTubuW7du\nqUGDBqnQ0FBVpEgR9cQTT6itW7eah69atUp5eXnZVa94eXllq3rFpEaNGubHKZs1a6a8vLyUwWCw\n+pgeG84pJ06cUG3btlX+/v4qNDRUDR8+XKWmppqHHzt2THl5eak1a9aY+x05ckR17txZlSlTRgUE\nBKiIiAi7R2iVUqpbt27q8ccfz3DZzZo1M6+nad0tt+XFixdVw4YNze14nnrqKbv2Cs5st/fee09V\nqVJF+fn5qeLFi6vGjRs7bMi7aNEi9cgjj6giRYqomjVr2jU2/eKLL1TFihWVr6+vCgsLUyNGjHD4\naPHbb7+twsLCHK7ztm3bVIMGDVRISIjy8/NTtWrVUmPHjrVqz7Fq1SrVqFEj8zhVq1ZVb7/9ttV+\n58y2sZSb1SuZ30YUDJFAXFxcXKaPoAmRl2zbts38iOn9uN8mJydToUIFxo8fn+1HKYUQ2ZPZ+cY0\nDIgCtt3rsvJC9YoQ4j4XHx/Pvn37qF+/PleuXOGDDz7AYDDQoUN2fopJCJHXSdIhhMgTxo0bx4ED\nB/D19aVu3bqsXbuW4sWLezosIYQbSdIhhPC4iIgItm7d6ukwhBA5LD++Bl0IIYQQ+ZAkHUIIIYTI\nFZJ0CCGEECJXSNIhhBBCiFwhDUmFyMP27dvn6RCEEAVcbp5nJOkQIg8y/X5Djx49PByJEOJ+YTrv\n5CRJOoTIg6pUqcLBgwe5du2ap0MRIsftOr+LY4nH+GzjZ0xuN5kaofa/fyNyVlBQEFWqVMl6xHsk\nSYcQeVRunACE8LQ0lUafyX2ILBvJ35/+TXCRYE+HJHKQJB1CCCFy3eHLhwF4qPhDLOq2iNIBpfHx\n9vFwVCKnydMrQgghco1Sil3nd/H8T8/zzu/vAFChaAVJOO4TUtIhhBAi18zcMZOXF73M0u5LqVeu\nnqfDEblMkg4hhBA57vTV01QoWoGYh2MoE1iGFg+28HRIwgOkekUIIUSOuXnnJnP3zKXKxCocTTxK\n4UKFafNQG0+HJTxEkg4hhBA5IuVuCvWm1mP3hd18/fTXPBjyoKdDEh4m1StCCCHcKjUtFYPBQOFC\nhRlYfyBNwprIuzcEICUdQggh3OhO6h2e/O5J/rPxPwC8XPdlSTiEmSQdQggh3GLh/oW0j21Pq8qt\nqFuurqfDEXmQVK8IIYRwizKBZSgVUIrBjw0mwDfA0+GIPEhKOoQQQrhs7Ym1dJvfjdS0VBpUaMB3\nnb6ThENkSJIOIYQQ2ZZyN4W3VrzFgUsHSLiRwNWUq54OSeQDknQIIYTItlSVyrIjy/D38Wd5j+UU\n8yvm6ZBEPiBtOoQQQjglTaXx5eYvebrK01QuXpm4l+Io5CWXEeE82VuEEEJk6fTV06w5voaJmydS\nyKsQA4oPkIRDZFte2GOaAMOBSKAs0An4OZPxnwVeBR4FCgN7gFHA8hyNUggh7mMT/5rI/H3z2dJ/\ni1SlCJflhTYd/sB24DVjt8pi/CeA34BodKKyCvgViMipAIUQ4n50NeUqfxz7A4BRzUbxV7+/JOEQ\n9yQvlHQsM36cNdSm+12gA9AeiHdXUEIIcT/bfnY7c/fM5dsd33J00FH8fPzw8/HzdFgin8sLJR33\nygsIAi55OhAhhCgITiSdoP60+oSHhLOhzwZJNoTb5IWSjnv1BhAAzPV0IEIIkZ8dTzpOWHAYYSFh\nLO2+lObhzfH28vZ0WKIAye8lHd2AEUBX4KKHYxFCiHxJKcX+i/upNakW38Z/C0DLSi0l4RBul59L\nOp4HpgKdgT+yGnnIkCGEhIRY9YuJiSEmJiZnohNCiHzi1cWvsjdhL9OfmU77au09HY7wkNjYWGJj\nY636JSUluXUZBrfO7d6lAR2BX7IYLwb4BngO/eRKZiKBuLi4OCIjI+89QiGEKCBup97G19uX9SfX\nk3grkXZV23k6JJHHbNu2jaioKIAoYNu9zi8vlHQEAFUsuiuhH3+9BJwCxgDlgJ7G4d2AmcAgYAtQ\nxtj/BiAv/xdCiCwopRixagQbT29k+QvLafxAY0+HJO4TeaFNRz109rQN/Y6O8cb/RxuHlwEqWozf\nHx33V8AZi8+EXIpXCCHyrYOXDhI1JYqKwRV5rtZzGPJcgbcoyPJCScdqMk9+ett0N8+5UIQQomAr\nF1SOKiWq8OSDT1K5eGVPhyPuM3mhpEMIIUQOOn31NNGzozmRdIJA30DmdJ4jCYfwiLxQ0iGEECKH\nzNg+gyspV0i5m8Llm5cJCwnzdEjiPiZJhxBCFGB//f0XPl4+/NEzyzcLiDysTx9Yvx4aN4bp0z0d\njesk6RBCiALmt8O/kXQriecefo6v2n4lL/nK5/r0gUWLICEBEhN1d35NPCTpEEKIAuLmnZvMiJ/B\nljNbuJZyjecefk4SjgJg/XqdcID+u369Z+O5F5J0CCFEAbHp9CbeWP4Gq3utpm65up4OR7hJ48a6\nhCMhAUJDdXd+JU+vCCFEPpaalsq8PfNQStH8weYcG3yM+uXr42WQ03tBMX06tGsHVavqv/m1agWk\npEMIIfKtE0kn2HBqA90XdGdj3408VuExSgeW9nRYIgfk50TDkiQdQgiRD6WpNJ7+4WnqlK3D/tf3\nU7VEVU+HJESWJOkQQoh8JPFmInfS7lAqoBSznp1FeEg4IUVCsp5QiDxAKv2EECKfuHTjEq1ntWbg\n0oEARJSJkIRD5CtS0iGEEPnAnN1zeGXxK0xuN5kGFRp4OhwhXCJJhxBC5GE379zEz8ePNg+14aPm\nH9G5Zmd5MkXkW7LnCiFEHpSm0lh6aCmVv6jM31f/JrhIMK/Vf00SDpGvyd4rhBB5TMrdFKJnR7M3\nYS996vShuF9xT4ckhFtI9YoQQuQxhQsVpk6ZOtQuXZunKj/l6XCEcBsp6RBCiDzgbtpdev/cmx92\n/QDA2JZjJeEQBY6UdAghhIdtPLWRWTtnUdi7sKdDESJHuSPpKAQ8DJwAEt0wPyGEuK9cvHGR+PPx\nLOm2hOAiwZ4OR4gc40r1yudAX+P/3sAaYBtwCmjupriEEKJAO3jpIO//8T5KKdpXa8/a3msl4RAF\nnitJR2dgp/H/9sCDQHXgP8BHbopLCCEKJKUU07ZNY8WRFczZM4dLNy8ByKOw4r7gyl5eAjhr/L8t\nMA84CMwAarspLiGEKJCupFxhxKoReBm82PXqLkr6l/R0SELkGlfadJwHagHngDbAAGN/fyDVTXEJ\nIUSBsuTQEuqVq0doQCi7Xt1FCf8Sng5JiFznSknHDGAOsBtQwEpj//rAPjfFJYQQBULy7WRWHVtF\nr4W9mLZtGoAkHOK+5UpJxyh0wvEAMBe4ZeyfBox1T1hCCJH/nbl2hk/Xf8pP+35iXZ91VClexdMh\nCeFRrj4y+5PxbxHAgC7x+NYdAQkhRH53N+0uTWY0YePpjSzptoRX6r5C1RJVPR2WEB7nSvWKNzAC\n+BtIRj+9AvAh6Y/SCiHEfen89fOMXDWSjac3UqRQEZqGN6VayWqeDkuIPMGVpONdoBfwFpBi0X8P\n0M8NMQkhRL6jlOJ40nEemvgQ5YuWZ13vddx89yb+Pv6eDk2IPMOVpKMn8BIwC+unVXYANVyYXxPg\nV3TJSRrQwYlpmqFfSHYLOGSMSQghPOL89fN4feBF/1/7M77VeHrU7kHjBxp7Oiwh8hxXko5ywOEM\n5uXjwvz8ge3Aa8ZulcX4DwKLgd+BR4EJwDSglQvLFkKIe3L66mnqTK4DQIWiFegf1Z+ihYt6OCoh\n8iZXGpLuA54Ajtv0/wc6eciuZcaPs14BjgDDjd0HgMeBocByF5YvhBAuGblqJHsv7mVg/YH0jOhJ\nuaByng5JiDzNlaRjNDATXeLhDTyLfg36i0A794WWoYakvxvEZDn6NexCCJHjvt/xPS8ufJF5XeZR\n3K84gx4bhMFg8HRYQuR5rlSv/Iz+zZWn0E+vfIBOOtqROyUNpdFvRbV0HigKyO9CCyFy1KQtk3hx\n4YsARJaNZHCDwZJwCOEkV39haC3QEghFt8l4HKnaEEIUYH9f/Zuu87pyJ/UO/9fw/0h+J5lKxSp5\nOiwh8hVXXw7mSeeAMjb9SgNXsX6E18qQIUMICQmx6hcTE0NMTIzbAxRCFCxd5nWhX51+HEk8wj8f\n/yeRZSM9HZIQbhcbG0tsbKxVv6SkJLcuw5UywbRMhil0Ow9XpQEdgV8yGWcs+tdtLX/R9gcgxNjf\nViQQFxcXR2SknCiEEM5LupVEz4U9+eXAL1QoWoGTQ05KVYq4r2zbto2oqCiAKPSrKu6JKyUdz9p0\n+wAR6HdljHJhfgGA5Q8SVDLO7xJwChiDbrRqehfH18DrwL/RPz7XAuiC44RDCCFcsub4Gg5fPsyu\n87sY32o8QxsO9XRIQuR7riQdCx30m4d+I+lz6HdmZEc94A/j/woYb/z/W6APuiqlosX4x4Gn0U+r\nDEYnJn2BFdlcrhBC2LmWco2F+xey8MBCDBg4OPAghbzyY020EHmPO4+kv4CpLky3mswbtPZ20G8N\nutpECCHcZs+FPTz834cBODX0FGUDy+LtdS81xkIIS+5KOvyBgehXmQshRL6SptJYeXQlX2/9mjKB\nZRjz5BgqFK3g6bCEKHBcSToSbboNQBBwA+hxzxEJIUQuOp50nP0X9xM9O5qVL6xk9rOz8fPx83RY\nQhRIriQdtq2p0oAEdPXK5XuOSAghcsl/t/yXAUsG8Fyt54h7KU4ehRUih7mSdHzr7iCEECI3XUu5\nxqrjqxiwZAAAo5uNplrJah6OSoiCz9mko3bWo5jtdCUQIYTIDalpqTSf2ZzapWuzrPsyWlVuJe/e\nECKXOJt0xDs53r2+HEwIIXLEndQ7FPt3MZLvJLO0+1KqFK9C5eKVPR2WEPcVZ5MO+YEBIUS+dSf1\nDm1/aEvynWQKeRXiqUpPyaOwQniAs0nH8ZwMQgghcsqSQ0sYsmwIwxsN5+MWH1O/fH1PhyTEfcvV\n93QYgBrAA4CvzbDMfjdFCCFyxb6EfdScVJOVL6ykbZW2dHukGwG+AZ4OS4j7mitJRyXgf8AjGQzP\n7O2iQgiR49adXMcTM54AYPeF3UxoM8HDEQkhwLUE4XN0dUspIBl4GGgCbAWauSswIYTIrtuptxmw\neAC7L+zmuVrPcWbYGQY3GOzpsIQQRq6UdDRE/7LrRfSLwVKBdcA/0QlJHbdFJ4QQTpq8dTInrpzg\nQvIF/H38+bHzj54OSQhhw5Wkwxu4bvz/Ivpn5w8AJ4HqbopLCCGckpqWyth1Y3lv1XsApI1Ik/du\nCJFHuVK9sof0l4VtBt4EGgMjgKNuiksIUcD16QPVqum/rjqWeIwJmybw58k/GVh/INffvi4JhxB5\nmCslHR+hf1UWdKLxK7AWuAQ876a4hBAFWJ8+sGgRJCRAYqLunj7d+emVUny64VOK+xXnyy1fEvdS\nHMX9iudcwEIIt3Al6Vhm8f8hdJVKCfSvz6a5IyghRMG2fr1OOED/Xb/e+Wkv3bjEw/99mHPXzzGp\n7SR2v7pbHoUVIp9wJenogX5kNtmi3yX3hCOEuB80bqxLOBISIDRUdztj1bFVbDu7jVt3bzGswTBe\nrfdqzgYqhHArV9p0TAAuAD8ATyO/tSKEyKbp06FdO6haVf/Nqmol6VYS+xL20XFOR1JVKgdeP8C4\n1uNyJ1ghhNu4UtJRFmgNdAPmADeBucBsYIP7QhNCFGTOtuHYcW4HEZMjANgzYA/VS1bHyyDvIBQi\nP3LlyL0DLEInHaWBIcCDwGrk6RUhhJukqTT2Juyl+czmAMx+djY1Q2tKwiFEPubqb6+YJAPLgWJA\nGPr3WIQQ4p7cTbvLu7+/y5w9c/jfc/+jQYUGFC5U2NNhCSHukau3DAHoBqVLgb+BocACoKab4hJC\n3KdeW/waPh/60OahNkzvMJ2m4U0l4RCigHClpONHoD26lGMu8AGw0Z1BCSHuP7dTbzNuwzgmbZ0E\nwIPFHiQ8JNyzQQkh3MqVko40oAv69eevIwmHEOIenbpyimpfVqOEfwm+7/Q9aSPSJOEQogBypaSj\nm92T1kUAACAASURBVNujEELcl27euYn/v/wZ8tgQutbsSqvKrSTZEKIAu9eGpEII4ZILyRdoPas1\nAKuOryL+lXgPRySEyGmSdAghct2UuCnsv7ifqLJRTGs/jahyUZ4OSQiRC+SBdyFErok/F0+dyXW4\nk3qHO6l3mNp+qiQcQtxH8krS8RpwHP12001AvSzGfxHYiX6C5gzwDSA/MSlEHvbz/p+pM7kO8efi\naVW5FRPbTpSfoRfiPnOvSUcRoKjNJ7ueA8YBI4E6wA7gNyA0g/GbAtOBKej3gnQB6gNTXVi2ECKH\nXb99nXd+f4cd53fwTLVnODzwMFVKVPF0WEIID3Al6QgAvgIS0CUNSRafRBfmNwydQMwE9gOvADeA\nPhmMXw9dKvIlcAJYb5y+vgvLFkLkoElbJrHt7Da+2f4NzcKbsfC5hVQuXtnTYQkhPMSVhqSfAC2A\nV4Hv0FUj5YGXgbezOS9fIBL42KKfAlYCDTOYZgUwCogGlgGl0KUdi7O5bCFEDrmbdpd/rvwn4zaO\nI8g3iPNvnMfPx8/TYQkhPMyVko72wADgJ+AusBb4CHiH7L/DoyTgDZy36X8BKJPBNDvQbTrmASnA\nWeAy+kVlQggPO3z5MD/t/YnJcZPpVL0TF4ZfkIRDCAG4VtJRHDhi/P8q6Q041wNfuyOoLDQAZqDb\ngPyGfjPqp8Zl98uF5QshHFBKseP8Dt5a+RYAxwYfo6R/SQ9HJYTIS1xJOo6if8r+JHAA3RB0M9AO\n3a4jOy4CqUBpm/6l0SUYjgxFJxvjjN270W1L1gLvYl9qAsCQIUMICQmx6hcTE0NMTEw2QxZC2EpI\nTqDUZ6UA2DNgDw8EP0Cgb6CHoxJCZEdsbCyxsbFW/ZKSsntZz5wrSce3QASwBhgDLEJXbfigG4Vm\nx20gDmgJ/GLs5wU8CXyRwTQGdKJiKc1imEMTJkwgMjIym+EJIbJy6NIhhv42FIDeEb2pGSo/Ni1E\nfuToRnzbtm1ERbnvXTquJB3jLf5fCVQHooBD6HdnuDK/mcBWYAswBPBDV6GATmzKAT2N3QvRic8r\nwHKgLDAB+As458LyhRAuuJ16m9XHVxM9O5qfuvzEtGemUSYwo6ZYQgjhWkPSF9Hv5zA5DsxHP+76\nogvzmwu8AXwAbAdqA23Qj+SCblBa0WL8H4DB6NKVXcbp9wHPurBsIYQLlh1eRuGPCvPniT+Z03kO\nHat3lIRDCJElV6tXlgG3bPoXNQ77zoV5fmX8ONLbQb//Gj9CiFx0O/U25caV49LNSwDULVeXjtU7\nejgqIUR+4c4ffCtP9huSCiHyCaUUhT8qbO5OeS8FX29fD0YkhMhvspN0bLf4fyXWjTm90U+0LHNH\nUEKIvEMpRbvYdiw5tITJ7SZzNeUqbzR6w9NhCSHyoewkHT8b/z6KfmQ12WLYbeAYum2HEKKAOJp4\nlMpfpL+2vHdEb3y8fTwYkRAiP8tO0jHK+Pc48CP2bTqEEAXI70d/p+X3LQEY2mAo41uPz2IKIYTI\nnCtPr3yLJBxCFFhXbl0heGwwhbz0Pcmmvpsk4RBCuIUrSUchYDj6nRrn0L8sa/pcdl9oQgh36tMH\nqlXTfzMSuyuWkH+HcDXlKsuPLEeNVDxW4bHcC1IIUaC58vTKCKA/+jXkH6J/ITYc6GjsFkLkMX36\nwKJFkJAAiYm6e/r09OFpKo2oKVHEn4sHYH2f9TSq2MhD0QohCipXko7u6KRjEfpH12LRPwC3E/1j\nbEKIPGb9ep1wgP67fn36sL9O/8W+i/vYdX4XAElvJRFcJNgDUQohCjpXko4ypL/u/DpgOjstRv/E\nvRAij2ncWJdwJCRAaKjuBhixagQf/qkLKG+9e4vChQpnMhchhLg3riQdp9G/hXIS/YuzrYFtQF0g\nxX2hCSHcZfp0XaWyfr1OOD76/AyG0eUBCC4czNaXtkrCIYTIca4kHQvRvwK7Cf1LsLOAvsADwH/c\nF5oQwp1MbTg2ntrIt/GrAIh+KJrF3RZjMGT4A81CCOE2riQdb1n8Pwdd4tEIOAj86o6ghBDul3w7\nmae+f4qNpzcyqe0kTgw5wQPBD3g6LCHEfcQdv72y0fgRQuRRC/cvpNOcTgDUDK1J38i+8rspQohc\n50rS8ST6Z+TDAEX668/XuDEuIYQbpKk0Jv41kSG/DQFgSrsp9I/q7+GohBD3q+wmHZPRj8smoqtT\nvIDGwOvAJONfIUQeEH8unjqT6wAQFhzG2t5rqRhc0cNRCSHuZ9l5I2knoBfQBygJNAQeA0KBfuhk\npIOb4xNCuOCVRa+YE46Yh2M4PuS4JBxCCI/LTklHb/TTKd/a9E8FpgPV0AnJzwghPOLGnRsE/CvA\n3L2p7yZ5jbkQIs/ITklHJPC/TIb/D/2uDiGEB7z060vmhKOkf0nSRqRJwiGEyFOyU9JREjiVyfDT\nQIl7C0cIkV237t7C72M/c/eW/luoW07yfyFE3pOdpOP/27vz+Kjqe//jr0mAFBGMIuIOggFpqyDU\nNbSlblj16nWrDa7kVysWi2hdsLRsrQtatBVUtBrQa03LrdoqFrDU9rZFBSWI1soqxA0lQFhKCEtm\nfn98zsycmcwkcyYzmZPJ+/l45EHONvP5zgxnPvmunYB9TRzf55wjIq1k0UeLGDpraGS77sd1dO7Y\nuYkrRERyx+volSlAXZJjXZLsF5EMC4VClEwvYVPdJgDuO/M+7hx6ZzNXiYjklpek4+9YZ9FkAmiu\nDpGsW7lpJcc9clxke/u47XQt6prDiEREUuMl6RiWrSBEJDUX/fYiXlr5EgDD+w5n/lXzcxyRiEjq\nMjENuohkWc3OGq5/+fpIwlH1/SpOPOzEHEclIuKNkg4Rn7t/0f3cudD6a1RcWMG1g66lIOBltLuI\niD8o6RDxqb0Ne/nS3V8iGAoCMPWsqYw8cWSOoxIRSZ+SDhEfevOTNxn5x5GRhGPHXTvYv9P+OY5K\nRKRllHSI+EgoFGLIE0NY9vkyAJ69+FmuPOHKHEclIpIZXpKOE1I879004hgN3A70BJYDPwTeauL8\nImACcCVwKLABm0NkVhrPLeILa7asoWR6SWR73c3r6F3cO3cBiYhkmJek450UzgkBhR5juAKYBtwA\nLAZuARZgc4LUJLlmDra6bTmwBjgsjecV8Y3Fnyzm1KdOBeCSAZfw/Heez3FEIiKZ5yXp6JOlGG4F\nngCedrZHAedjCcXUBOefC3wDOAbY6uz7KEuxiWTVtvptFE8tBuArPb7Cg8Mf5Jy+5+Q4KhGR7PCS\ndKzPwvN3wlavvdu1LwQsBE5Lcs2FwNvAOOAqYCfwEvBToD4LMYpkxZz353DF76+IbL8z6h06FKib\nlYjkLy+D/UuA3wIHJDhWDFQCxyU41pSDsWaRL+L2b8T6aiTSBxgKfBn4b2AscBnwqMfnFsmZi393\ncSThePS8RwlNDCnhEJG85+Uudzu2tP22BMe2OsfGAde1PKwmFQBBrBPpDmffrcDvgRuB3YkuGjt2\nLMXFxTH7ysrKKCsry16kInFeXfsqw58dHhn+WntnLcVfKm7mKhGR7KusrKSysjJm39atW5OcnZ6A\nh3NXYc0ZS5Ic/xrwHNDPw2N2wppHLsWaSMKeBroBFye45mngdKzmJWwA8L6zb23c+YOBpUuXLmXw\n4MEeQhPJrCtfuJLn3nsOgHvOuIe7vn5XjiMSEWlaVVUVQ4YMARgCVLX08bzUdBxF42YQt03OOV7s\nAZYCZxFNOgqAM4GHk1zzT6w5pQuWsIAlOkHgE4/PL5J17s6iACtvWkm/7l5ycxGR/OClT8c2oG8T\nx/sC29OI4UHgeuAarMbiMaAz0Tk37iU6sgWsNmWzc3wANpLlAeApkjStiOTKmHljIgnHkd2OJDgh\nqIRDRNotLzUd/wDGAK8lOT7GOcer8JwbU7DOo8uwYbHhOToOJbYGZSdwNjAdG8WyGfgd8JM0nlsk\nK+r21tHlni6R7Zrbazh4v4NzGJGISO55STruAd4EngfuB1Y4+wcAd2CJwulpxvGI85NIohWuVgKa\nzEB8KX4o7M4f72S/jvvlMCIREX/wknQswzp8zqJxB89NwOVY/wyRdikYCnLZnMt4ccWLAMw8fyY3\nfO2GHEclIuIfXicGmAv0AoZjI0UCWK3Dq0BdZkMTaTve/eJdBs4cCMBZfc7ihe+8QNeirjmOSkTE\nX9KZjagOeDHTgYi0VRf/7mL+sOIPAJzT9xwWXLUgxxGJiPiTl9ErZwAfYPNnxCsG/o316xBpFzbu\n3MgNL98QSTj+deO/lHCIiDTBS9IxFluYLdGw2K3ATOCmTAQl4nf3/fM+ev6iJ09UPcEDZz9AcEKQ\nrxzylVyHJSLia16aVwYBdzZx/M/YVOkieatubx09HuhB3V7rwvT4BY/z/SHfz3FUIiJtg5ek4xBg\nbxPH92HzbYjkpQ07NnD4g4dHtut+XEfnjp1zGJGISNviJen4FDgeWJPk+PHAhhZHJOIzoVCILvd0\nYde+XfTv3p+rT7ia8d8Yn+uwRETaHC9Jx5+wWUPnA7viju3nHJubobhEfGHV5lX0n9E/sv3ata9x\neNfDm7hCRESS8dKR9G7gIGxejjuBi5yfcc6+g5xzRPLCix+8GEk4LhlwCaGJISUcIiIt4KWm43Og\nFHgUmxI94OwPAQuA0c45Im3all1b6H5/98j2a9e8xreO+VYOIxIRyQ9eJwdbD5yH1WociyUeq4Et\nmQ1LJDd+8fovuP3PNgjryG5H8tHYjwgEAs1cJSIiqUhnRlKwJGNJJgMRyaWGYAP73bMfexr2ADD2\nlLE8dO5DOY5KRCS/pJt0iOSNyvcqGfHCiMj29nHbtW6KiEgWKOmQdi0wOdp08sqIVziv5LwcRiMi\nkt+UdEi7VL21mt6/6h3ZXvy9xZx8xMm5C0hEpB1Q0iHtzvi/jOc37/0GgHOPPZd5V87LcUQiIu2D\nkg5pNz7b8RlHPHhEZHvdzevoXdw7dwGJiLQzXiYHE8mZ8nLo39/+Tcf0xdNjEo768fVKOEREWpmS\nDvG98nKYOxdWrbJ/m0s83AlKQ7CB+xfdz5j5YwA45YhTCE0MUdShqBUiFxERNzWviO8tWgQ1NfZ7\nTY1tJxNOUGpq4KOv/pBZP5sBwB2n38G4oeM4sPOBrRCxiIgkoqRDfK+0FGprLZHo0cO2k4kkKHd1\no75oBwBf7vFlpp49tXWCFRGRpNS8Im3C7t3QrRtccAFUVCQ/r+Rbb8CkADgJB/N+yZ6H3m+dIEVE\npEmq6RBfCzeXbN9utRxNGfDIAFYctiK6Y1IQCFDdMashiohIipR0iK+l0p9j++7tHHDfAZHtbqtG\nsaPyMUJAIAC9erVOrCIi0jQ1r4ivlZZGazgS9ecY/5fxMQnHxts2su03j9G3L3TsCH37wurVrRiw\niIgkpZoOaTUlJVBdbTUPqSYCFRXWxLJokSUc4f4cwVCQwimFkfMO6nwQm+/YHNluSaKR6PlERKTl\n/FLTMRpYD+wC3gROSvG6UmAfsCw7YUmmlJTA2rWwd6/9W1KS+rUVFbByZTQBmPn2zJiEY/mo5TEJ\nR0t4nRNERERS54eajiuAacANwGLgFmAB0B+oaeK6YuAZYCFwSJZjlBaqroZQyH4PhWw7He5VYQFC\nE0MtjCyWlzlBRETEGz/UdNwKPAE8DawARgF1QHN/Y84EngXeAALNnCs51quXdeqE9Dp3flj7YUzC\nMf/K+RlPOKD5PiQiIpK+XCcdnYDBWG1FWMjZPq2J60YCvYHJKOFIqKVrlWTa6tWk3bkzMDlA34f7\nRrYbJjQw/NjhWYjSmnAuuAD69Wt+ThAREfEm180rBwOFwBdx+zcCxyW5pgS4FxgKBLMXWtvlngq8\ntta2/fDl6bVzZ93eOrrc0yWyPfqk0cw4b0aGo2rMD6+ViEg+ynVNh1eFwHPARGBNjmPxrXzol/Cj\nBT+KSThqbq9JKeHwWw2PiIhE5bqmYxPQAPSM298T2JDg/K7AEGAQEP4GKsCaWPYCZwN/S/REY8eO\npbi4OGZfWVkZZWVlaYbuX17WKvGbUChEwZRoLnzmMWey8JqFTVwR5dcaHhGRtqCyspLKysqYfVu3\nbs3oc/ihP8SbwBJgjLNdAHwEPAzcH3duABgQt280cAZwKTbsti7u+GBg6dKlSxk8eHDmova5tjjX\nxIsfvMglcy6JbL876l2O73l8ytf3729DXcP69bOhtiIikp6qqiqGDBkC9gd/VUsfL9c1HQAPYiNX\n3gbeAsYCnYFZzvF7gcOBa7FOpv+Ou74GqE+wv11rK4lGWCaGwrblGh4RkfbAD3065gC3AVOwSb5O\nAM4lOkfHocBRTVwfcn6kDVqxaUVMwjH929PTHgqrkSciIv7mh5oOgEecn0RGNnPtZOdH2pjL//dy\nfv/v30e2d/9kN50KO7XoMZVoiIj4l1+SDmlH4ofCntP3HBZctSCHEYmISGtQ0iGtavY7sxn5x2jl\n1bqb19G7uHfuAhIRkVajpENaTbbXTREREX/zQ0dSyXNvffpWTMJReWmlEg4RkXZINR2SVYNmDmL5\nF8sj28EJQQKB5qeHaYvzjIiISNNU0yFZsWXXFgKTA5GE46aTbiI0MZRywjF3rk30NXeupjQXEckX\nqumQjCutKOX1j1+PbF9TXc/0iUUpX58Pa8eIiEhjqumQjNkX3EdgciCScBRs6wOTQsybW+SptqK0\n1GYUBc0sKiKST5R0SEZULKug4886RrZ7v7SG4ENrAe+1FZpZVEQkP6l5RVos0VDY8mqY+3H666Bk\nI9FQ51QRkdxSTYek7ds3/i0m4Vh49cLIUFi/1Vaoc6qISO6ppkPSEpgcsKX4HNdVN3Bmn9gcNteJ\nhps6p4qI5J5qOsSTmp01sc0pb42CSSFeX+Tvj5I6p4qI5J5qOiRlF//uYv6w4g+R7e5PbmbzJwe1\niS/xigr16RARyTUlHdKsYChI4ZTCyHZRYRH1P6mnvLptfYm3hRhFRPKZkg5p0ph5Y5i+ZHpke1H5\nIk4/6nRAX+IiIuKNkg5JSqvCiohIJvm795/kxJJPl8QkHFOGTVHCISIiLaaaDokRX7ux96d76VCg\nj4mIiLScajoEgJ17dsYkHAN7DiQ0MaSEQ0REMkZJh3DXwrvY/979I9urf7iad0a9k8OIREQkH+nP\n2Dzgnn8CYoexNjU3RSgUomBKbN6Z674bmktDRCR/Kelo48JritTUQHU1BAJQXw+1tVBSAtu22bE1\na+Af/4DVq+26P63+E+c/d37kcV696lXO7nt2jkph3GWprbVtJR4iIvlDSUcb515TZPfu6P6aGti6\nFfbute1gENatsy/yWb38ORRW66OIiOQ39elo40pLoaio8f4ePaBXLyhwvcMN3d+LSTgePe9R3yQc\noPVRRETynZKONq6ionHS0bGjLSe/ejX06QOFhcBtPeEHJ0TOqR9fz40n3di6wTajosLi7tfP/lXT\niohIflHzSh649NJoX4gePWK/sN/7oJ7Od3eOnNu9c3c23bEpR5E2T4mGiEj+Uk1HHkhWQzBm3piY\nhKN6bLWvEw4REclvfqrpGA3cDvQElgM/BN5Kcu4lwI3AQKAIeB+YBLya9Sh9Kr6GQOumiIiI3/il\npuMKYBowETgRSzoWAD2SnP915/i3gcHAX4GXgUFZj9Tn5rw/JybhOGnLNPo9F6K8PIdBiYiI4J+a\njluBJ4Cnne1RwPlAOTA1wfm3xG2PBy4C/gtot1NpxtduXFcd5JW5Ac17ISIivuCHmo5OWG3FQte+\nkLN9WoqPUQB0BTZnNrS24e3P3o5JOPoc2IfQxBCvLwrEzHvx7LPQqZNNGiYiItLa/FDTcTBQCHwR\nt38jcFyKj3Eb0AWYk8G42oShFUNZ9HF0Fq3t47bTtagrYPNc1NZawlFQEJ0obO1aSzzCs5OKiIi0\nBj8kHS01ApgAXAi0m6EZDcEGOvws9u2L7yzqXnvFnWCEQjZluoiISGvyQ9KxCWjARq249QQ2NHPt\nd4FfA5cBrzV14tixYykuLo7ZV1ZWRllZmadg/eDSOZfywgcvRLaXj1rOCT1PSHhuuA9HSYnVcIRC\ntj5Lr16tEamIiLQVlZWVVFZWxuzbunVrRp8j0PwpreJNYAkwxtkuAD4CHgbuT3JNGfAUNvLl5SYe\nezCwdOnSpQwePDgz0eaQu+9G991D2HTP2ylfW1JiNRy9eqlpRUREmldVVcWQIUMAhgBVLX08P3Qk\nBXgQuB64BhgAPAZ0BmY5x+8lOrIFrEnlGeBH2Fwehzo/3Vop3lb3xsdvxI5O+eQUCp58m5IS6N+f\nlIbErl4Ne/Yo4RARkdzwQ/MKWAfQHsAULHlYBpwLOGMvOBQ4ynX+9VjC9IjzEzYbG2bbZoX7YJSW\nRptG+k3vx+otrkxhcgOECqgBNm+2FWQ1JFZERPzOL0kHNE4g3EbGbX8ry7HkRHl5dA2V2loY8b1N\nVB4VnR9t9EmjWXDzDNZiY4rBEg7QUvAiIuJ/fmleESxpiMyrMeiOmITj41s+ZsZ5MygosM6gYeGl\n67UUvIiI+J2fajravdJS2FIbYtPo2Fyw33MhvvWsJRjBoCUY4RVlDzjA9rubY0RERPxISUcOJBtF\n0umSUWzq9Xhk+8LPqpj75ImsCkbPCScaBx6oRENERNoWJR2tzD1fhntm0ESrwnbqFO2zEVZTYwnH\nypWtGLSIiEgGqE9HK6uujvbJCIVg/ZaPYxKOh4Y/FJlZdN++xter74aIiLRVSjpaWcxMoJMC7Btz\ndGSz60N7mTh8bGTOjcLC2Gs7doQLLlCTioiItE1KOnKhwy6YFK3d6L3zcgqmhNixrQPbt0NlJXTr\nBg0Nrks62MReSjhERKStUp+OVrZmYBlc9dvojseWs2XXCTF9N+rr7cct4JcJ60VERNKkmo4MKi9P\nPiX5yPKQ9d043kk4PjkZJoXgi8YLtQUC0LVrNNHQAm0iIpIPlHRkSHg20VWr7F934nHZ99cwu5fr\npZ79V3hyMWBzbBxyCBQVRQ/37Qvbt9u/HTvav1ovRURE2jo1r2RIzGyirinJA5MDcITrxEmhmOuC\nQUs8RoxovOaKEg0REcknSjoypLTU1ksJzxTa94y/E5j8zegJs/8K64c1ui48BFYdREVEJN8p6ciQ\nioroCrGrywqZF3D1DL1nB+zZv9E1GgIrIiLtiZKODJr56z0U/bwodmdcc0pY167Wb0NERKS9UEfS\nDOl7x4jYhONXHyZMODp2hJEjlXCIiEj7o5qOFiovh1m9AtDFtTNJ7QbAMceoOUVERNon1XS0wJXf\nq2XWga6hKW/c0mTCoXVTRESkPVNNR4rCnUTDI00Om3YYnx/1uR1cNhL++BSQfNpQdRoVEZH2TklH\nCsITf9XUwJpN1cya3Dt68H/mw9rhja7p2hV69rRVZXv10pwbIiIial5JIH468+efdyb+OvEpgmN6\nR0+8b0ujhGPkSFuyfvt2SzT27FHCISIiAqrpaMRdq1FbCyUlsH1HECa51pmvPwDu2xrZ7NgRrrpK\nTSciIiJNUdIRJ3468xreh4lfjZ5w/0ao6xHZLCqyKcyVcIiIiDRNzStxSkttlAlAoOsXMNpJOD48\nw0amuBIO0JLzIiIiqVLSEaeiArp2/w+Bb/6cUN2BBP40A6Z9As/8JeH59fXRxd1EREQkOSUdCew7\nYDWhU6fBoe8QWjIadhzR6JzwUvSae0NERCQ1SjoSOPPLJ9L9f6rh05Pp0cOGv7ode6z14+jXT3Nv\niIiIpEodSROwFWO7sWgRBJ3FYnfsgMJCm8ZcQ2BFRES8U9KRRLj2on//6GiWhgYoUN2QiIhIWvzy\nFToaWA/sAt4ETmrm/GFAFVAPrAauzVZg7tEsfum/UVlZmesQMiafygIqj5/lU1lA5fGzfCpLpvkh\n6bgCmAZMBE4ElgMLgB5Jzj8GeAX4CzAQ+CXwJHBONoKrqLB+G37qv5FPH+h8KguoPH6WT2UBlcfP\n8qksmeaH5pVbgSeAp53tUcD5QDkwNcH5o4C1wO3O9kpgKHAL8Go2AvRDoiEiItLW5bqmoxMwGFjo\n2hdytk9Lcs1pceeDJRvJzhcREREfyHXScTBQCHwRt38jcGiSa3omOP8LoBtQlNHoREREJGP80LzS\nKj744INch5AxW7dupaqqKtdhZEQ+lQVUHj/Lp7KAyuNn+VSWTH935nrlkE7ATuBS4CXX/qexmouL\nE1zzf9jIlVtc+0YCDwHFCc4/DHgLaDytqIiIiDTnU2xU6YaWPlCuazr2AEuBs4gmHQXAmcDDSa55\nAzgvbt/ZwOtJzt+AvViHtShSERGR9mkDGUg4/OI72Pwc1wADgMeBzUSHzN5LdGQLQG/gP9jIluOA\nHwB7scRDREREpEnhycHqsZoM9+Rgs4DX4s7/JrGTg12T/RBFRERERERERERERERERERExCPfLhiX\nJi/luQT4Mzah2jZsFE9W1qFJk9f3JqwU2Acsy05YafNaniLgbqJ9ltZhQ7z9wGtZrgHexYa5fwY8\nBRyUxfhS9Q3gZWxYXxC4KIVrhuHfe4DX8vj9HpDO+xPmt/tAOmXx8z0gnfL49T7Qaq7A3shrsdEs\njwNbaHrBuJ3AA0B/7Ma7F//8J/VanoeA24AhQF/sw70bGJT1SJvntSxhxdj6OvOxLwa/SKc8f8S+\nBM4AjgZOAU7Pbpgp8VqWb2I3/5uAXtiXwXvA81mPtHnnAlOA/8ZunBc2c77f7wFey+PnewB4L0+Y\nH+8D6ZTFr/cA8F4eP98HWs1iYuf0CACfAHcmOX8qlqW5VQLzMh9aWryWJ5F/AT/NZFBpSrcsvwUm\nYysP++UvHPBennOBWhJPWpdrXstyG7Ambt8PgY8zH1qLpHLj9Ps9wM3Ll7SbX+4B8byUx6/3gbBU\nyuLne0C8VMrT4vtArtdeaal8WzAunfLEKwC6YnOd5FK6ZRmJzcUymdzPmOuWTnkuBN4GxmFfBbv7\niwAACUFJREFU6Cuxv66/lL0wU5JOWf6MrYf0bex96QlcDrySvTCzxs/3gEzwyz2gJfx6H/DKr/eA\ndLX4PpDrGUlbqqkF445Lck1zC8btzmSAHqVTnni3AV2AORmMKx3plKUEmwxuKJZ1+0k65emDlWUX\nVn3ZA3gU6A6UZyfMlKRTluVYW+7/YklLB2wW4ZuyFGM2+fkekAl+uQeky8/3Aa/8eg9IV4vvA229\npkNijQAmYLO8bspxLF4VAs9hVanx1XdtVQF207wS+2tnHnAr1o+ira2IfCo2Ud9ErJbkXKxvxMxc\nBiWNtOV7AOTffSCf7gGQgftAW6/p2AQ0YH+5uPUk+Tzxn2PVQ/Hnbyf3f+GkU56w7wK/Bi6j8Qyu\nueC1LF2xjnCDgBnOvgKsCi88zf3fshFoitJ5bzZgvbt3uPatwMp0JNZJLhfSKcstwAJgmrP9L6wz\n5j+A8TSuOfAzP98DWsJv94B0+P0+4JVf7wHpavF9oK3XdLgXjAsLLxj3RpJr3nCOuzW1YFxrSqc8\nAGVABXbT8UtnOK9l2QZ8FRjo+pmJtYEOBJZkM9gUpPPe/BM4HKvqDuuH/eXzSRZiTFU6ZQlgiYpb\n0HWsLfHzPSBdfrwHpMPv9wGv/HoPSFc+3QfSlm8Lxnktzwgs/huxv97CP91aKd6meC1LvEn4q9e6\n1/J0AT7C2tYHYGPiVznX5Vo6n7M9wCisnboUeIumk+HW0gX7y3gQdgMc6/x+lHO8rd0DvJbHz/cA\n8F6eeJPwz33Aa1n8fA+A9D5rfr0PtKp8WzDOS3n+imWewbifitYINAVe3xu3ifhnfH6Y1/L0x0ZG\n7MRuPg/gn7Zcr2W5kWh16qfAM8BhWY+yecOIfu7d/xfC/wfa2j1gGN7K4/d7wDC8vz9ufroPDMN7\nWfx8DxiG9/L49T4gIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIl4U5joAkXZm\nHRACFuc6EBEREfGX2cCLGXy87kBn13YQuDCDjy/ty2wy+/kUyaq2vrS9SLaFnJ9M2ZxgX1tdnbED\nsC/XQWRRIZYUZvL9FxERSWo2sX9J/g34FXA/lkBswBakcpsEVGOLiX3qnB+2Hhjj+t29QNeHSZ4T\n4JfYwl7JTKLxSpxjseacsGHY0uD/AWqxZbePdh2/CFtYaxewFphAbBNsEFtd8iXnMSYkiWU9MB5b\nCGqHE8P5QE/gZWffcmBI3HVDgX8AddjCWL8C9nMdvxp4G9iOve6/IboqLsCBzr6NzmOsAq5zlT1I\n7Mqr4ZU1w6/Bddjr8l/Av7GVW4/GFuf6BbYU+X+AN7EF44i77nxsCfadwO+wGq1yp/xbnPIUuK5L\n9XHPAT7AXrd52AqyYO95/CJv3wA6ATOAz7D3cj0wDhER8b3ZNE46tgI/BfpiX4QNwFnO8cuc48OB\nI7GVW/+f6/p1RJOOg7EvimuAQ7CmF7CVHV+Ii+OXNL0S5ySaTjo6OHFNBY7BVr68mugS1l93jl+N\nLf1+FpYEuROLIPA5cK1zzlEkth7YBFyPvUaPYF+erwKXAiVO+d53XdMX+1Id4/x+GrCU2JVSR2Kv\na2/gFGAR8Irr+AwsaRqMJQtnAhc4x4aRWtKxG0t8TnXi7Az82tlXir12P8K+zI+Nu24+MBB7LWuc\n8v4WOA5LSOqB77ieP9XHXeCU6UTnNXvWOd7FefxXsM/PIUBH4DYs6S3F3qNS4ApERMT3ZtM46fi/\nuHMWA/c6v98KrCB506U76YDEfTrinxNaXtNxENG/hBNZCNwZt+8qrKbGHeu0JmIIWwc87dru6Vw7\nybXvFGffIc72k8DMuMcZijXfdEryPF9zHiNcG/JH4Kkk5w4jtaQjCBzvOudorMYjfunuPwN3x113\njOv4Y1jthbumZp6zvyWPeyNWyxM2m8aflV9h76eI76hPh4g3IeC9uH0biH55zgFuxmoJ5gN/wpoU\nGlorwCS2YF9QC7AvtoVYrJ87xwcCpwM/cV1TiDUBfAn7Kx2seSMV77p+3+j8+16CfYc4vw/Evuyv\ndJ0TcH6OwZothmCJywlYU0oB9n4cjSV6jwHPY7UCrwJ/AN5IMd6wPXFxHo+9DqvizivCanPC6oht\nytrobNfF7Qt/TtJ93M9dj5HMbOw9Xol9Buc62yI5p6RDxLu9cdshop1BP8GaLs4CzgYeBW7H2upT\n7XQZpHHn0o4ZuKYceBg4F6tu/7kT42Ksqn4CjZt1wKr4w3Y2E0eY+zUKNbEv3MehC1bT8XCCx/rY\nOb4Aqy0YgTVf9HL2hWtC5jv7zsPK9Resaed27PWB2Nco0Wu6K257fyxhHEzjxPE/rt8TfSbi32/3\n56Qljxv/Psd3dF2GJWrfxj6Hc7Ak83JEckxJh0jm1WN/Xc7FvvRWAF8F3klw7l4az5ezEfhK3L5B\nxH75x6sh2sHQfU38F9I7zs99wOtAGZZ0VGF9Dz4kN6qwMid7/hOwJqJxRJt8Tk5w3iasA+szWH+J\nB7Cko8Y5fjiwzfl9UApxLcPen55Yx9tMydTj7iHxfXwHlmzMAX6PJWTFWL8dkZwpaP4UEXEJV/kn\n23cdVqPwVaAP1jGzDuvYl8h67K/RQ7EmA7AOo19zri0BJtM4CYn3V2wkxx1YR8zRWI1G2DFYv5NT\nsdqAc5zH/sA5PgXr0DrBea4BwHeBnzXzvJkyFWvemY4lAyXYaJrpzvGPsC/YMdjreiHWmddtirP/\nWKwM4VEoAGuwGpNJzvHzsY6bzVmFjYh5BrgYex1PBu7CalS8Cn9OMvW467CErB/WMbkD1q/ou1gS\n2Q/rvLoBJRziA0o6RJoWP09Honk73PtqsVEb/8SGhZ6BffnVJnn8H2FNAR9hozXA+iP8DBuWuwRr\nWnimmThXAD/Ako13sKTlF67jO7Fmn+extv7HsdEeT7ie8wIsGVmC9YW4GUuKssX9Or6HNUH1A/6O\n1XxMJlqrUYMldJdjIzjuwF4792PsxhKr5Vhn373Yly/O72XYF/G7WO3HeBK/l/FGYq//NOx1fhHr\nX1LdxHXNfU7Sfdz4fb/G3s+3gS+wkSrbsdfnLey9PJr0EiQRERERERERERERERERERERERERERER\nERERERERERERERERERERERERERERERERERERERERERERaS3/H6W/OZpU6vMLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x48ac490>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# lists to hold data\n",
    "mine = []\n",
    "cci = []\n",
    "\n",
    "with open('output_cleaned_cci_db_4.csv', 'r') as datafile:\n",
    "    csvFile = csv.reader(datafile, delimiter=\",\")\n",
    "    _ = next(csvFile)\n",
    "    for row in csvFile:\n",
    "        mine.append(float(row[3]))\n",
    "        cci.append(float(row[4]))\n",
    "\n",
    "\n",
    "mine_np = np.array(mine)\n",
    "\n",
    "cci_np = np.array(cci)\n",
    "cci_np[cci_np == 9.96921e+36] = np.nan\n",
    "\n",
    "coeefs = np.polyfit(mine_np,cci_np, 1)\n",
    "#r squared\n",
    "R2 = str(trendline(mine_np,cci_np, Rval=True))\n",
    "print(R2)\n",
    "m, b = np.polyfit(mine_np,cci_np, 1)\n",
    "\n",
    "plt.plot(mine_np,cci_np, '.', label=\"matchup values\")\n",
    "plt.plot(mine_np, m*mine_np + b, ':', label=\"regression R2 = {}\".format(R2))\n",
    "plt.ylabel(\"CCI Data values\")\n",
    "plt.xlabel(\"Insitu user measurements\")\n",
    "plt.title(\"Regression of CCI daily values against Insitu User Values\")\n",
    "plt.legend(prop={'size':10})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by/4.0/\"><img style=\"float: right\" alt=\"Creative Commons Lizenzvertrag\" style=\"border-width:0\" src=\"https://i.creativecommons.org/l/by/4.0/88x31.png\" /></a>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}