GEE_histograms.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "______\n",
    "<img style=\"float: right;\" src=\"./images/headerlogo.png\">\n",
    "## Google Earth Engine Tutorial\n",
    "### Histograms\n",
    "_____"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In previous examples, multiband images were displayed in RGB combination providing the *min* and *max*  parameters to strecth digital levels to 0 - 256 values.  A Styled Layer Descriptor (SLD) can also be used in GEE for image display. <br> \n",
    "With the RasterSymbolizer element we can specify the pixel strecth we want to apply. <br>\n",
    "There is only three possibilities:\n",
    "<ul>\n",
    "<li>Normalize: stretch the image contrast by mapping the minimum and maximum image values to the minimun and maximum brightness.\n",
    "<li>Histogram: stretch the image contrast by distributing an equal number of pixels at each brghtness value\n",
    "<li>GammaValue: adjust the brightness of the image. A value of 1 means no change \n",
    "</ul><br>\n",
    "More complicated visualizations can be created. For additional info, check the link: http://docs.geoserver.org/latest/en/user/styling/sld/reference/rastersymbolizer.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"https://earthengine.googleapis.com/api/thumb?thumbid=0cd951a04a68411bca99167f5e51313a&token=b2840a4fc3c3bdad9e71bf87f4b1dfb0\"/>"
      ],
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import ee\n",
    "from IPython.display import Image\n",
    "ee.Initialize()\n",
    "\n",
    "image = ee.Image(\"users/rosamaguilar/tutorial/subset\")\n",
    "\n",
    "# Define a RasterSymbolizer element with '_enhance_' for a placeholder.\n",
    "template_sld = '<RasterSymbolizer>' + \\\n",
    "    '<ContrastEnhancement><_enhance_/></ContrastEnhancement>' + '<ChannelSelection>' +\\\n",
    "      '<RedChannel>' + '<SourceChannelName>b5</SourceChannelName>' + '</RedChannel>' + \\\n",
    "    '<GreenChannel>' + '<SourceChannelName>b3</SourceChannelName>' + '</GreenChannel>' + \\\n",
    "    '<BlueChannel>' + '<SourceChannelName>b2</SourceChannelName>' + '</BlueChannel>' +\\\n",
    "    '</ChannelSelection>' + '</RasterSymbolizer>'\n",
    "\n",
    "# Get SLDs with different enhancements.\n",
    "equalize_sld = template_sld.replace('_enhance_', 'Histogram');  # Equalized\n",
    "normalize_sld = template_sld.replace('_enhance_', 'Normalize');\n",
    "\n",
    "# Display the image with a histogram equalization stretch.\n",
    "# band names of the outputs images are related to the type of sld used\n",
    "Image(url = image.sldStyle(equalize_sld).getThumbUrl({}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"https://earthengine.googleapis.com/api/thumb?thumbid=5c367efb1f49d19a93840d4bf44b4fc9&token=5628a6cc874157fb7d897c92b355c92c\"/>"
      ],
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Display the image using  Normalize stretch\n",
    "image = image.sldStyle(normalize_sld)\n",
    "Image(url = image.getThumbUrl({}))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Above examples shows a predefined style, but it could be of interest to modify the histogram applying a customized stretch. First we should explore the original histogram. To do this, we will apply a *reducer* called *frequencyHistogram* over a region that for this case will be the whole image. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A **reducer** makes common calculation such as minimum, maximum, median, first, count, etc. It can be applied to a neighbourhood, region, etc. \n",
    "\n",
    "The function *reduceRegion* will be applied over the whole image to obtain the values for the histogram. \n",
    "\n",
    "The following parameters can be specified for the function <i>reduceRegion</i>:\n",
    "\n",
    "+ *reducer*: the *reducer* to apply.\n",
    "+ *geometry*: region over which the reducer will be applied. By default is the footprint (extent) of the image's first band. <br> \n",
    "+ *scale*: scale in meters of the projection. \n",
    "+ *crs*: the projection to work in.\n",
    "+ *crsTransform*: a list of values to apply a transformation of the CRS. Default value: null.\n",
    "+ *bestEffort*: a flag to indicate that the calculation can be done at larger scale to succeed if too many pixels are involved. Default value: False.\n",
    "+ *maxPixels*: the maximum number of pixels to process. Default value:10000000.\n",
    "+ *tileScale*: A scaling factor to aggregate tile size. Default value:1. Using higher values would reduce the memory needed to finish the computation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"https://earthengine.googleapis.com/api/thumb?thumbid=27475a6f49961e81335e287011c043a1&token=061cab463b162d083fa3cfb8d17a2355\"/>"
      ],
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# visualize the image without any stretch\n",
    "# load the image\n",
    "image = ee.Image(\"users/rosamaguilar/tutorial/subset\"); \n",
    "# Visualize the image\n",
    "Image(url=image.getThumbUrl({'min': 0, 'max': 2048,'bands': 'b5,b3,b2'}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Minimum and maximum value for : b5 is 0.00 9114.00\n",
      "Minimum and maximum value for : b2 is 0.00 10486.00\n",
      "Minimum and maximum value for : b4 is 0.00 8015.00\n",
      "Minimum and maximum value for : b8 is 0.00 10267.00\n",
      "Minimum and maximum value for : b3 is 0.00 8531.00\n",
      "Minimum and maximum value for : b1 is 0.00 7679.00\n",
      "Minimum and maximum value for : b7 is 317.00 8766.00\n",
      "Minimum and maximum value for : b6 is 10.00 9437.00\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "# compute the histogram \n",
    "# set the geometry of the image to use as region in the function reduceRegion\n",
    "geom = image.geometry()\n",
    "# set the original scale of the image\n",
    "scale = image.projection().nominalScale()\n",
    "# compute the histogram\n",
    "bfreq = image.reduceRegion(\n",
    "  reducer = ee.Reducer.frequencyHistogram(),    #.unweighted(),\n",
    "  geometry = geom,\n",
    "  scale = scale,\n",
    "  maxPixels = 10000000 ).getInfo()\n",
    "#print(bfreq)\n",
    "\n",
    "# The output of reduceRegion is a dictionary with the format:\n",
    "# {'bandname': {'key':value}}\n",
    "    \n",
    "# We will build series  for x values and y values in order to build the histogram    \n",
    "\n",
    "# define two dictionaries for frequencies and values respectively\n",
    "freqs= dict()\n",
    "vals = dict()\n",
    "\n",
    "# populate dictionaries\n",
    "for l in bfreq.keys():\n",
    "    \n",
    "    d = bfreq[l]\n",
    "    #print(d)\n",
    "    freq = []\n",
    "    val = []\n",
    "    for k,v in d.items():\n",
    "        #print(k,v)\n",
    "        if not(k == 'null'):\n",
    "            val.append(int(k))\n",
    "            freq.append(int(v))\n",
    "    freqs[l] = freq \n",
    "    vals[l] = val\n",
    "    print('Minimum and maximum value for : %s is %.2f %.2f' % (l, np.min(val), np.max(val) ))\n",
    "    #print(freqs[l]),np.min(val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhAAAAF7CAYAAACD27KNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3X24HHV99/H31wS1VoyioKgEfLhAtKWYiBat1lKt3mot\ngk9HRIui+FgbwKp4ewtafKpNan0oScWCEE8tF8HaWlArEQVBIQcSREiCEBJkIydCQiAEA/nef8wc\ns+dwnuZkz87u2ffruvba3fnNzH53J9n9nN/M/CYyE0mSpCoeUncBkiSp+xggJElSZQYISZJUmQFC\nkiRVZoCQJEmVGSAkSVJlBghJklSZAUKSJFVmgJAkSZUZIKQOFhEfiYilddfRzSLizoh40RhtH4+I\nC9pYywciYnm7Xk+aTgYIqSYRsTwi/maU6Tsj4hCAzPx0Zh4ziXV9PCKWTUed7RQR6yJiW0TcFRG3\nR8R/RcTTpvll2z2ev9cP0IxggJA6T8f9wETErDa9VAJvyMxHAQcAW4Cvtem1JVVggJA6T/zuwYgu\n9oj4bEQ0ImJLRNwQEa+IiL8CTgFeFRFbI+Kuct7ZEfHpiLglIn4dEf0R8bimdT0zIi4v1/WDct3L\nm9p3RsR7I+Ja4O6I+P2IWBARa8oegrUR8d6m+fcvlzkuIn5ZzvOZiHhCRHyvfJ3lEbHPZN5/Zm4D\nvgHMG9YYMS8iLo6I35S1HN/UFhHxyYjYGBG3RsR7JvF57xERXy3rWx0RRzat76URcWVEbI6IX0XE\nlyPiYU3tN0fEB8vP8a7y/T1prM8YeOIk6pG6ggFC6nwJxY8Z0AccmplzgJcAazLzP4FPAf+dmXuW\nf71DESpeATwfeEo5bWm5rtnAt4HvAI8FPgK8jQf3fvSVr/MoYBuwDnhx+RrHA/8QEYePWObFwLOA\n5wJ/C3wT+BvgccCOsq4JRcQc4FhgddO0xwPfA76cmY8FjgROjYg/K2c5DngL8ELg6cBzgEdO8FIv\nB64AHgOcBPRHxNDntQ04PjMfDbygfG8njlj+GOAN5fvbBnyyrHUWxWf8fYrP+KMUn5k0IxggpHp9\nJiLuaLrdydi7MHYADwP+MCJmZ+atmXnjOOt+M/DJzPxV+df8icBLI+IJwOHAXsCnMvP+zPwZxQ/9\nSJ/NzF9n5o4sXJCZtwFk5iXAdyl+VJt9MjO3Z+YNwErg0sy8ITN3ABcwokdhFEvLz+FO4FCKEDHk\nWOCSzDy/rOEXwFnAm8r2NwFfzMy1mbkd+DAw0e6X1Zn51czcmZn/DSynCE5k5mWZubJ8vA5YMsr7\n/XJmrs/M31IEtPnl9OdTBIfTys/4Ckb/jKWuZICQ6vXhzNyr6fYYmnZhNMvMHwIfp/gLdzAizouI\n/cdZ95OBW5qWbwD3ldP3BRqZubNp/vWjrGND85OIOCYiVpS7D+4E/g/FX96/exng9qbn24Bfj3g+\nUY/Am8rP4SBgNtB8EOUBwCtHBK73A08o25/I8Pd8O8V7Hs8tozx/EkBEHBYR3y93iWwGTmf4+2XE\n+7sH2LN8vC9wW2Y+MM5rSV3LACF1kcw8IzMPB+YCvwX+uWzaOcrst1L84AJQ9jw8tJx+G/CEiGj+\nDpg7yjp+t96I2I/ir/2TgceVP/IXMkbg2Q1Dx0CsBT4InNF03MEGYFlz4MrMOZn5l2X7bcDvQlV5\nvMXDGN/IEDYX+FX5+BvAxcAB5W6MjzL593sb8MQRB6CO9hlLXckAIXWJiHhORBweEXtQ/FV9D3B/\n2fxrYP8RP1bnAqdExJMj4pHAQuD7mbmRYp//ZuAj5cGWhwGvn6CEoZ6DwbKeVwB/MbLMKb69UWXm\nt4BNwPvKSecAR0TEUWXde0TEH0XEc8r2fuC9EXFgRPwexbEhDzx4zcMcGBFvj4hZEfFK4M+Afy/b\n9gQ2Z+b2iDgYeHeF8q8A7gA+Vtb5PIpjJaQZwQAh1WesYx3Gmv4o4CsUP6i3UXSRf6BsOw+4i2LX\nxh3ltE9THKNwOXATxbEAxwJk5v3AXwF/SfEj9xmKH+fm7v5hdWTm9RRd+MvLGl4H/OcEtVc9JXW0\n+T8DfDAifq88/uJlwAlAo7x9iV27Db5GEZx+DNwIDABbJ3jNi4A/pvgcFgHHZOZNZdsJ5WvfRfHZ\n90+i3qKh+IxfTXGQ5m8owsyZE9QidY3I7LhTziXVICLOoPhOOKHuWiR1PnsgpB4VEX9S7t6IiPhz\nijMY/qPuuiR1h1oCREQ8NCK+WA4CszIivl5Of3pEXFYO5vLTiHhm0zIj2w6uo3ZpBnkqxX76rcCX\ngQ9l5g/qLUlSt6hlF0ZELAIekpkfKJ/vk5m3lyO1nZWZ50TE0RRfaM8t5xmzTZIktVfbA0REPILi\nwKcnZebdTdP3BtYCew2dmx4RDYrR37aO1dZ0sJMkSWqTOnZhPI3iaOePlmPMXxIRRwD7MfrANnMn\naHuQiHhEFOPlP2J63oIkSTPTZH9DZ7eroBGvuT/w88z8SEQcSjFW/Ctp3TnkhwKXAQMRcfeItoso\nTm2TJKnXvYziVONmj6QYcv4FwE/GWrCOALGeYmCXbwBk5jURcTNFqHhCRDykqadhbjn/VmDfMdpG\nc0B5P9qY+y+iOB9bkiSN7QA6KUBk5m/KAyJfDlxYXvXuAOBSikFfjgXOjojXAhuGjnGIiDHbRrEO\n4Nxzz+Xggz1Zo5MsWLCARYsW1V2GRnC7dCa3S2ea6dvl+uuv581vfjOUv6VjqaMHAorhYM+MiM9S\n9Ea8MzMbEfEu4KyIOAXYQnFp3iHjtY20HeDggw9m3ryJLvyndpozZ47bpAO5XTqT26Uz9dB22T5e\nYy0BIjNvBo4YZfoaikvgjrbMmG2SJKm9HIlSkiRVZoCQJEmVGSDUVn19fXWXoFG4XTqT26UzuV0K\nM/JqnBExD1ixYsWKXjnQRZKklhgYGGD+/PkA8zNzYKz57IGQJEmVGSAkSVJlBghJklSZAUKSJFVm\ngJAkSZUZICRJUmUGCEmSVJkBQpIkVWaAkCRJlRkgJElSZQYISZJUmQFCkiRVZoCQJEmVGSAkSVJl\nBghJklSZAUKSJFVmgJAkSZUZICRJUmUGCEmSVJkBQpIkVWaAkCRJlRkgJElSZQYISZJUmQFCkiRV\nZoCQJEmVGSAkSVJlBghJklSZAWKyGg049dTiXpKkHmeAmKxVq+C004p7SZJ6nAFisjZtGn4vSVIP\nM0BIkqTKDBCSJKkyA4QkSarMACFJkiozQEiSpMoMEJIkqTIDRFXLljmYlCSp5xkgqjJASJJkgJAk\nSdXVEiAiYl1EXB8RV0fEQES8rpz+9Ii4LCJWR8RPI+KZTcuMbDu4jtolSVJ9PRA7gddn5rMzc15m\nnldOXwyckZkHAZ8DzmpaZmTb2e0s+EG8uJYkqYfVFSCivO2aELE3MB9YCpCZ5wP7RcRTx2tra9XN\nGo3i4loGCElSD6rzGIivR8TKiPjXiHgcsB/QyMydTfOsB+ZO0CZJktpsdk2v+8LMvDUiZgGnU+yO\n+BgjeiV214IFC5gzZ86waX19ffT19bXyZSRJ6kr9/f309/cPm7Zly5ZJLVtLgMjMW8v7ByLin4DV\nFD0K+0bEQ5p6GuaW07eO0zamRYsWMW/evNYUfeedrVmPJEkdYrQ/qgcGBpg/f/6Ey7Z9F0ZEPCIi\nmrsF3gQMZOYmYAVwbDnfa4ENmXlTZg4CA6O1ta3wzZvb9lKSJHW6OnogHg+cHxEPodhlcRPw1rLt\nXcBZEXEKsAU4rmm58dra64Yb4JJLant5SZLq1vYAkZk3A6PuV8jMNcDzq7a13cqVsGRJ3VVIklQb\nR6KciosvrrsCSZJqZYCYiquuqrsCSZJqZYCQJEmVGSAkSVJlBojdNThYdwWSJLWdAWJ3bdpUdwWS\nJLWdAUKSJFVW17UwukejAYsXw7331l2JJEkdwx6IiQxdtnvr1rorkSSpYxggJElSZQaI3bVsWdFL\nIUlSDzFA7C4DhCSpBxkgJElSZQYISZJUmQFCkiRVZoCYrCuvrLsCSZI6hgFisryEtyRJv2OAkCRJ\nlRkgJElSZQYISZJUmQGiVRoNOPVUB5WSJPUEA0SrDF10ywAhSeoBBghJklSZAUKSJFVmgJAkSZUZ\nICRJUmUGCEmSVJkBQpIkVWaAqIvjRkiSupgBoi6OGyFJ6mIGiE5gb4QkqcsYICYyODj9r2FvhCSp\nyxggJrJpU90VSJLUcQwQrdCOXgpJkjqIAaIVlizZFSIaDVi//sHzrF8/+vTBQVi8eHrrkySpxQwQ\nrbBs2a5dHUcdBQcdNDwsrF9fTBs5HYrllixpX62SJLWAAaJV7ryzuP/tb2H79uHHTmzaVEwbOV2S\npC5lgGiVzZvrrkCSpLYxQLSKAUKS1ENm113AjLFt2/DnCxfC4x8PJ59cTz2SJE0jA8R0Wbq0uD/m\nmHrrkCRpGrgLYyxDw0uvW9e+11y82NEoJUldodYAERHHRcTOiHh1+XzviLgwItZExKqIeGHTvGO2\nTYuh4aV/9avpWfdoYz8sWWKAkCR1hdoCRETsDxwPXN40+TPA5Zl5IPA24BsRMWsSbfW78srJz9to\nOPaDJKmr1RIgIiKArwLvA37b1PR64AyAzLwKuA340zHaftXUVr+rrhp9eqMxfq/CsmXTU48kSdOo\nroMoTwR+nJlXF1kCImIvYHZm3t403zpg7hhttwBzp73Se+7ZveVf8xrYuXPsdgOEJKkLtT1ARMSz\ngKOB6T2GoVV+85vdW37HjuHPBwfHDw1emEuS1AXq6IF4IbA/sLbclfEEYAlwKnB/ROzT1NNwAHBL\nZt4REaO1jXJ1ql0WLFjAnDlzhk3r6+ujr69v8tVu3z75eSdj7drxj39wqGtJUpv09/fT398/bNqW\nLVsmtWzbA0RmnkF5LANARCwHFmbmf0XEc4F3A6dFxGHAE4EflbOeN0rbJeO91qJFi5g3b97uFfzL\nX+7e8iNt2NDa9UmSNEWj/VE9MDDA/PnzJ1y2EwaSSiDKxx8GzomINcB9wDGZ+cAk2qZPq8eB2Lix\nteuTJKkGtQeIzDyi6fHtwMvGmG/Mtq6yu8dUSJLUARyJst0mOqZi2TIHk5IkdTwDRLvdd9/47QYI\nSVIXMEC02913112BJEm7zQDRbjfcUHcFkiTtNgNEu7V6XAlJkmpggJAkSZUZICRJUmUGCEmSVJkB\nohN5QS1JUoczQHSixYsdC0KS1NEMEJ3oggvg1FMNEZKkjmWA6FRLlhggJEkdywAhSZIqM0BIkqTK\nDBCSJKkyA4QkSarMADEWx2KQJGlMBoixbNpUdwWSJHUsA0QnsxdEktShDBCdzF4QSVKHMkBIkqTK\nDBCSJKkyA4QkSarMANHJbrqp7gokSRqVAaKTrVpVdwWSJI3KANHJ7ryz7gokSRqVAaKT3Xdf3RVI\nkjQqA0Qnu+MOWL++7iokSXoQA0Qn+8Uv4KCDDBGSpI5jgOh027c7IqUkqeMYICRJUmUGCEmSVJkB\nQpIkVWaAkCRJlVUOEBHxrclMUwsNDtZdgSRJw0ylB2LuKNOeuruFaByehSFJ6jCzJztjRJwAvAs4\nMCIGmprmANe1ujBJktS5Jh0ggIuA1cC/AAuapt8FeNWn6eQ1MSRJHWbSASIzbwFuAQ6evnI6RKMB\nS5fWXcUuGzbUXYEkScNU6YEAICIOAD4EPK15+cw8omVV1a3RgAsvrLuKXTZurLsCSZKGqRwggP8A\nfgB8CXigteVoVNu21V2BJEnDTCVAPDwzP9LySjS2e++tuwJJkoaZymmcP4+I0U7l1HS57rpit4ok\nSR1iKj0QewMrI+JyYPvQxMw8arIriIjvAo8HkuIsjg9k5jUR8XTgbOBxwGbguMz8RbnMyLa/zszr\np1B/91m3Dlatgn33rbsSSZKAqQWIc8vb7nhdZt4FEBFHAmcBhwKLgTMy85yIOLqc/txymZFtZze1\nzXwOJiVJ6iCVA0Rmnr27LzoUHkqPBnZGxN7Ac4CXlvOcHxFfioinAluB+aO1ZeZNu1tPV3AsCElS\nB5nKaZxfG216Zr6t4nrOBv6MYjfGK4D9gNsyc2fTbOsphs6+C2iM0dYbAWLZMjj6aHdjSJI6wlR2\nYaxoevxw4GhgYIx5x5SZbwWIiGOBzwEfA2IK9YxpwYIFzJkzZ9i0vr4++vr6Wvky7bF8eXEgpQFC\nktQi/f399Pf3D5u2ZcuWSS0bmblbLx4RjwDOy8xX7sY6tgH7AzcCjxnqaYiIBvACil0Ya4G9RraN\ntgsjIuYBK1asWMG8efOqFzQwAPPnT/HdTKMVK2Aq70eSpEkaGBhgfvEbOD8zx+wgmMppnCNtB54y\n2ZkjYk5E7Nv0/EhgU2YOUvRuHFtOfy2wITNvKtsGRmtrQf3dw8t6S5I6xFSOgVjY9HQWxYGPv6iw\nijnAeRHxcIrjH24HXlW2vQs4KyJOAbYAxzUtN15bb/BMDElSh5jKMRDNO0fuB/4ZWDbZhTNzPfC8\nMdrWAM+v2iZJktprKqdxnjYdhUiSpO5R+RiIiNgzIr4cEWvK25ciYs/pKE6SJHWmqRxE+RWKnovX\nA68r1/GVVhZVq0YDFi6ceL46OJiUJKlDTOUYiEMy84+anr8nIla2qqDaNRqwdGndVYxuw4a6K5Ak\nCZhaD8Ss5l0W5eNZrStJY9q6te4KJEkCptYDcTZwRUR8s3z+euDfWleSxnTPPXVXIEkSUCFARMSj\nKEaC/IeI+Dnw52XTV9j9q3NqMn70I4ezliR1hCq7MD5HcUVMMvPCzDw5M0+mGGb6s9NRnEZYt64I\nEJIk1axKgHhRZp4/cmJmngO8qHUlaVwLFxoiJEm1qxIgHhinbec4bWqlpUsNEJKk2lUJELPL4yCG\niYg5wB6tK6lm3XDBKgOEJKlmVQLEvwPnRMRjhiaUj/+tbJsZuuGCVa95DaxfX3cVkqQeViVA/D2w\nGdgQEVdHxNXABoqDKD85HcVpDDt2dEfQkSTNWJM+jTMzHwDeGhGfAOaVkwcy85fTUpkkSepYU7ka\n5y8BQ4MkST1sKkNZS5KkHmeAkCRJlRkgRvKS2ZIkTcgAMdLmzXVXMDndMF6FJGnGMkB0qxUr6q5A\nktTDDBDdavVqOPVUR6WUJNXCANGttm2D004zQEiSamGA6FbXXFN3BZKkHmaA6FY33lh3BZKkHmaA\nGKlbzsKQJKlGBoiRtm2ru4JqPJ1TklQDA0S386qckqQaGCAkSVJlBoiR7rmn7gokSep4BoiRuu0Y\niPPPdywISVLbGSC63QUXGCAkSW1ngBjp3nvrrkCSpI5ngBhp+/a6K5AkqeMZIGYCx4KQJLWZAWIm\nWLu27gokST3GADET/OxndVcgSeoxBoiZ4Oab665AktRjDBAzwX331V2BJKnHGCAkSVJlBghJklSZ\nAWKkbtwdsGNH3RVIknqMAWKkbgwQd91VdwWSpB5jgJgJuu0CYJKkrtf2ABERD4uICyLihoi4OiK+\nGxFPK9v2jogLI2JNRKyKiBc2LTdmW0t14+6AzZvhyivrrkKS1EPq6oFYnJnPyMxnA98GvlpO/yxw\neWYeCLwN+EZEzCrbPjNOW+vcc0/LVznttm+HF7wA1q+vuxJJUo9oe4DIzPsy86KmSVcA+5ePXwec\nUc53FXAb8Kdl2+tHtP2qqa117r+/5atsix074NJL665CktQjOuEYiA8A34qIvYDZmXl7U9s6YO4Y\nbbcAc9tXZhdYubLuCiRJPWJ2nS8eEacATwPeCTyi1etfsGABc+bMGTatr6+Pvr6+Vr9UZ9i4se4K\nJEldpL+/n/7+/mHTtmzZMqllawsQEXEycCTw55m5HdgeEfdHxD5NPQ0HALdk5h1jtI2703/RokXM\nmzdvmt5BB7rttrorkCR1kdH+qB4YGGD+/PkTLlvLLoyIOBF4I/DSzNza1HQe8O5ynsOAJwI/Gqft\nkpYW1mjAnXe2dJVtNcnUKEnS7qrjNM4nAZ8H5gDLy1M5Ly+bPww8PyLWAF8DjsnMBybR1hqNRnFK\nZLdavx5OOql4H5IkTaO278LIzF8xRnApd0+8rGqbSr/+NSxcCMccA/vuW3c1kqQZrBPOwpAkSV3G\nADETLV7sbgxJ0rQyQMxES5YYICRJ08oAIUmSKjNANBscrLuC1plJ70WS1HEMEM02baq7gtaZSe9F\nktRxDBCSJKkyA4QkSarMANGsm4exHmnlymJkSkmSpoEBolk3D2M90sKFcNBBhghJ0rQwQDSbSWMn\nPPAAbN8Op58+s96XJKkjGCCa3X77xPN0GweVkiRNAwNEs3vvrbsCSZK6ggGi2fbtdVcgSVJXMEA0\nu+++uiuQJKkrGCCGNBpw0011VyFJUlcwQAxpNOC22+quQpKkrmCA6AULF3omhiSppQwQvWDpUgOE\nJKmlDBBDVq+uuwJJkrqGAWLIL39ZdwWSJHUNA4QkSarMACFJkiozQAzZurXuCiRJ6hoGiCEbN9Zd\ngSRJXcMAMWTbtrormF6Dg3VXIEmaQQwQQzZvrruC6XX55XVXIEmaQQwQQ7ZsqbuC6fWJT8C8eXDN\nNXVXIkmaAQwQvSITrr4arruu7kokSTOAAWLIjh11V9AeDmstSWoBAwQUP6i33lp3Fe1x4YUGCEnS\nbjNAQPGDumlT3VW0j2dkSJJ2kwGiF/VSWJIkTQsDRC869VTPxpAk7RYDRC+68Ua49NJi182pp3pM\nhCSpMgNEr9q8uQgOp50Gq1YZJCRJlRggVBwTcdppBghJ0qQZIKA3z0rw6qOSpN1ggIDePCth40ZY\nvLjuKiRJXWp23QWoJv/zP70ZnCRJLWEPRK8yPEiSdoMBQnDnncX9okXFgZSe3ilJmkDbA0REfCEi\nbo6InRFxSNP0p0fEZRGxOiJ+GhHPHKft4JYV1GjA+ee3bHVdacOG4v7cc3cFCM/KkCSNo44eiPOA\nFwDrRkxfDJyRmQcBnwPOGqft7JZV02jABRe0bHVd6eKL665AktRl2h4gMvPSzLwNiKFpEbE3MB9Y\nWs5zPrBfRDx1vLZ21z5jXXXVrseLF/fmaa2SpEo65RiI/YBGZu5smrYemDtBm1ptyRIPsJQkTahT\nAoQ6ydBBlZIkjaFTxoHYAOwbEQ9p6mmYS9HTsHWctnEtWLCAOXPmDJvW19dHX19f6yqfiTZvLu4X\nLy7Oxth331rLkSRNj/7+fvr7+4dN27Jly6SW7YgAkZmDETEAHAucHRGvBTZk5k0A47WNZ9GiRcyb\nN286S5/ZliyBE04wQEjSDDXaH9UDAwPMnz9/wmXrOI3zjIjYADwJ+G5ErCmb3gWcEBGrgb8Djmta\nbLw2SZLUZm3vgcjMd40xfQ3w/Kptmgaf/eyux4ODxW4MeyIkSU08iFIPdvfdux57qW9J0igMEJIk\nqTIDhH9Zj2/olE4Hl5IkNTFADJ2yqNENXSfjfe+Da66ptxZJUscwQGh8GzcW9zfeCNddV28tkqSO\nYYDQ+JYt2/X4zDPhpJPc7SNJMkBoAs1nZCxfDgsXGiAkSQYISZJUnQHCC0dV13xGRqNRDDRlr4Qk\n9RQDhGdhVLd27a7Q0Gg40JQk9aCOuJhWrbZurbuC7nPiibBjBxx+OOy9d93VSJJq0Ns9EI0GXHRR\n3VV0nx07ivtNm+qtQ5JUm94OEKtWFTdNjcePSFLP6u1dGGvX1l1Bd/vHf4THPrbuKiRJNejtHggP\noNw969YV40JAcX/SScVw156VIUkzXm/3QGj3XXVVcb90aXE/b15xVsarXw377ltfXZKkadXbPRCS\nJGlKDBBqraFrZwwOuitDkmaw3g4Q115bdwUzz1CAWLLEAaYkaQbr7QCxZk3dFcxczVfxBIe8lqQZ\nprcPonQcg/a45hp4+9thYGBXkPAAS0nqar3dA3HPPXVX0Buuu64ID1Ds2rAXQpK6Xu8GiEYD7r67\n7ipmvuuum3ieoV6J0caQcNeHJHWk3g4Q27fXXcXM99a3wkc/Onza4GAx6NRJJw2/oud11xX3q1Z5\ntU9J6nC9GyAGB+uuoDdkwi23DJ+2dm0xcuXChUVYWLy4mD504OWmTbsfGuy5kKRp1bsBwitJ1mfD\nhl2P164tjouAB5+58Y53wA03jL+u9euL20j2XEjStOrdALFuXd0V9K6vfGXX409+8sHtP/pRcT8w\nACtXFo9H6zFavx4OOqi4jRYiJEnTpncDxNAPk9qv+eDV229/cPtQjwTARRcV90M9Rs27JjZtKo5j\n2b4drr9+eLu7qCRpWvVugLBruzusWlXcL106/kGVa9cWwWHVqqLdXVSSNK16M0A0GhPvW1dnufBC\neN/7hvcsND/evLkIDmvXFs8dJEySplVvBogf/tC/ULvRsmVw/PHF49NPhy98YVfbhRfumgeKQAHF\nGR5DPReelSFJLdObAeLSS+uuQFN1663F/bJlu0IDwE9+UtwvX17cD7UtWQLz5hXze1aGJLVMbwaI\n0Q7c08wyFCgANm6c3IiYkqRJ680Ace+9dVegdhu67olnZ0hSS/Te1TgbDf8a7UU331zcf+EL8L3v\nwcknF88///nifrTnXjFUksbUewFi1SoHkepFQ2dnXHhhcbv99mKY7aVLi+l33w1HHVUMrw1wzDEG\nCEkaR+8FiK9+te4KVIdf/3r483PPHf58yRL4wz9sXz2S1OV66xiIRqM4hVO9J3PieT71qV2Ph07/\nBE8BlaRRzOwAMfKAuVWrHP9BY2sOCEuW7BoF83//tzgF9IMfNERIUmlmB4jmsNBoDB94SJrIpk3F\n8RJveUvxfOlSeyIkqTSzA0SzVauGDzwkTeSv/xre857h05YscWRLSaJXAsQ118AJJ9RdhbrN/feP\nfsbO6afDy142fGRLA4WkHjOzA8SyZfD978MrXgG33FJ3NQL66y6gFZYtg2uvLR6///1wyCFwxBG7\nAkUXhon+/hmxZWYct0tncrsUuipARMTTI+KyiFgdET+NiIPHXWDZsq77Ip/pZtx/u5/8pAgTQ1d3\nveIK+Iu/KMLEq161a9yJF7+46AnrUH4hdia3S2dyuxS6bRyIxcAZmXlORBwNnA08d9wlmq+JIE23\nD3yg2PVvJ2NcAAAID0lEQVQBMDBQ9H4NufRSOPtseM5z4F//FT70IXjc4+A73yl2sTlwlaQu0jUB\nIiL2BuYDLwXIzPMj4ksR8dTMvKne6qTSUHgYzfvfP/z5JZfserx4MSxaBFdeCXfdBe98Jxx22PTU\n2Arr1xf3c+d25/oljW2S1wzqmgAB7Ac0MnNn07T1wFzAAKHutnEj9PXtej40YuqeexZB4he/gJe8\npDil9MMfhm3bYI894PDD4fjj4fGPL55/7nPwpjcVyzQaRTAZ6t0Y+lHeY49d1/x44xuL+9Wri/lH\n9oIMreOVryxeY+7cYj0HHbRrOShG+hzqSYHhrzsZzYHhyiuL9zV7NqxZ8+AQYbiQptckx0vqpgBR\nxcMBrq+7Cj3IFmCg7iK6ydatcPHFxeNzz33wENzXXlucWgrwyEcW1/T4p3+COXNgy5Zi+umnF6Fh\ntKvQltf+2AIMPOMZxTruuw/e8Y4iEOy5Z7Eb8BOfgFmzimDy85/D9u3F8s98ZrHezOJ29tnF8o1G\nUdesWcV7eMxj4ClPKS5qtnlzERBe9Sr44heLdd16Kzz0oXDGGcXtgQeK20EHwd/+LdxxR3Gtkvvv\nhze8AX77WzjzTNhnn+JYp6OOgr33LgLN5z8Pb387rFxZTIfh84xm9Wr49KeLXUqbNsFHPgJ77bVr\nucmsYxps2bKFgYFx/scMvd+TT94V6iZjcBDOOad4fOyxw9/T4GD19zqVZXZHu19vhAm3SxU1v5fR\nXD908cHyt3QskZMZ4rcDlLsw1gJ7DfVCREQDeMHIXRgR8SZgafurlCRpxjgmM78xVmPX9EBk5mBE\nDADHAmdHxGuBDWMc//Bd4BhgHbC9fVVKktT1Hg4cQPFbOqau6YEAiIgDgbOAx1L0uh6XmdfVWpQk\nST2oqwKEJEnqDF01kJQkSeoMMy5AVB6tUlMSEQ+LiAsi4oaIuDoivhsRTyvb9o6ICyNiTUSsiogX\nNi03pTZVFxHHRcTOiHh1+dztUqOIeGhEfLH8HFdGxNfL6SO/s57ZtEyVNr/rpiAiXhERK8rvsVUR\n8ZZyuv9fJpKZM+oG/AA4tnx8NPCzumuaiTfgYcDLm56/F1hePv4a8P/Kx88BNgCzyudnTqXNW+Xt\nsz9wWXl79e589m6Xlm2TRcAXmp7vU96P+Z011TZvlbbLb4BnlY/3B+4Fft/vsUl8dnUX0OJ/CHsD\nm4GHNE1rAE+tu7aZfqMYJfSm8vHWoS/H8vlPgSOm0HbFUJu3StsigO8DzwaWNwUIt0t92+QRFAd+\nP3LE9DG/s6baVvd77bYbMAj8Sfn4kPIHfw//v0x865rTOCfJ0Srr8wHgWxGxFzA7M29valsHzJ1C\n2y0U207VnAj8ODOvjggA3C61expwB/DRiHgJsA04jSIEjPWdddcU2/yuq+aNwAURcQ/waOAoYE/8\n/zKhGXcMhNovIk6h+II8pe5ael1EPIuiO/v0umvRMLMpusd/npmHUQTub5bTo87CellEzAL+L3Bk\nZh4AvAQ4F7fLpMy0ALEB2Dcimt/XXIpkrmkQEScDR1IcD7E9M+8A7o+IfZpmOwC4ZYptbrtqXkjx\nQ7U2Im4G/hhYArwet0ud1gMPAN8AyMxrgJspttUTxvjOGu/7zO+61jgU2DczLwPIzKuAWyl2Zezw\n/8v4ZlSAyMxBikstHAswwWiV2k0RcSJF999LM3NrU9N5wLvLeQ4Dngj8aIptTZes1EQy84zMfFJm\nPjUzn0Kx//UdmXkGbpfaZOZvKA56fDlARDyF4oflUsb4zhrv+8zvupYZCmLPgOLMFopjTG7A/y8T\nq/sgjFbfgAOBnwCrgZ9RHl3rreWf85OAnRTXJxkArgYuL9v2oRgCdQ1wLfCipuWm1OZtytvpYnYd\nROl2qXdbPKXcHqvK/y9HltPH/M6aapu3StvlDU3bZCXwhnK6/18muDkSpSRJqmxG7cKQJEntYYCQ\nJEmVGSAkSVJlBghJklSZAUKSJFVmgJAkSZUZICRJUmUGCEmSVJkBQtKkRcR3IuI9o0y/JiKOHGe5\nj0fEwumtTlI7GSAkVXEm8LbmCRHxHOAJwH/VUpGkWhggJFXxbeDJEfEHTdOOA74OHBwRP46IqyLi\n5+Vl3h8kIt4aERc0PX9lRCxvev7miLiiXM8PI+KQcvrzymkDEbEqIk6YnrcoaTJm112ApO6RmfdH\nxDnA24EFEfEwoA84nOIyyEdk5o6IeDjwk4j438z82WirGu15RDy/XN8Ly/X8CcUlsP8A+AjwD5n5\nzXLeOdPwFiVNkgFCUlVfA34YER8EjgZ+kZmrI2Jv4F8i4lCKK7U+GTiU4kqRk/VXwCHATyMiymmP\nLoPKcuBjEXEgcHFmXtai9yNpCtyFIamSzLweuBF4NcXui6+WTZ8CBoE/ysxDgUuAh4+yivuBWU3P\nm+cJ4OzMnJeZzy5vT87M+zLzC8ArgduAT0XEl1v6xiRVYoCQNBVfA04BDgP+o5z2GODWzMyIOAh4\n6RjL3ggcEhEPi4jZwJua2r4NvDki9gOIwvzy8YGZeUtmnkkRVp7X8ncladLchSFpKr4JLAL+PTO3\nldP+HjgnIt4IrAN+MNqCmfnTiPgf4FrgZmCAMgxk5qUR8XfABRExC3go8B1gBfC+iDgC+C1FL8aJ\n0/TeJE1CZI48lkmSJGl87sKQJEmVGSAkSVJlBghJklSZAUKSJFVmgJAkSZUZICRJUmUGCEmSVJkB\nQpIkVWaAkCRJlRkgJElSZQYISZJU2f8HFcltjEQIPhkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1cbc73652b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhcAAAF7CAYAAABhB6n0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xm8XXV97//XR6KiglE0ShwCimLRNtJE2+LVlmIdqq31\ngXUISHv1+pM69No49Dr0/gra24p6Sa1DSa5SkSFVrwFtEdEqUJlJDpCASQhD5g3nnAyHM0/53j/2\nStwJZ9xZe6+1d17PxyMP9l7fNXy+ex9y3vmu71orUkpIkiTl5XFFFyBJktqL4UKSJOXKcCFJknJl\nuJAkSbkyXEiSpFwZLiRJUq4MF5IkKVeGC0mSlCvDhSRJypXhQmpREfGpiLi86DpaWUTsiYjfnaTt\nbyPiyibW8pGIuK5Zx5MayXAhlVBEXBcR/32C5fsiYiFASukfUkpnz2BffxsRqxpRZzNFxOaIGIiI\nRyOiMyL+LSJOavBhm/18BJ/HoLZguJBaS+l++UTEUU06VALemVJ6KnAi0ANc3KRjS5oFw4XUWuLA\ni0OG7SPigoioRERPRGyIiDdFxJ8Anwb+KCJ6I+LRbN05EfEPEbElIh6JiJUR8cyafb00Im7J9vWz\nbN/X1bTvi4gPRcQ6oC8inhIRSyPivmxkYVNEfKhm/ROybd4TEQ9k63w+Io6PiJ9kx7kuIp41k/6n\nlAaAK4BFBzVGLIqIn0fErqyW99W0RUR8LiIejojtEfHBGXzej4+Ib2T1bYyIt9bs73URcUdE7I2I\nHRHxtYh4Yk37QxHxiexzfDTr33Mn+4yB58ygHqklGC6k1pag+osOWAKcmlKaC/wBcF9K6QfA3wP/\nnlI6NvtXP1QDx5uAVwEvyJZdnu1rDvBD4GrgGcCngPfy2FGTJdlxngoMAJuB07NjvA/4YkScdsg2\npwMvA34L+CvgO8B/B54JjGZ1TSsi5gLnABtrlj0b+AnwtZTSM4C3AudFxO9nq7wH+DPgNcCLgFcA\nx0xzqDcCtwJPBz4GrIyI/Z/XAPC+lNLTgP+S9e2jh2x/NvDOrH8DwOeyWo+i+hn/lOpn/Bmqn5nU\nFgwXUnl9PiJ21/zZw+SnRUaBJwK/ERFzUkrbU0r3T7HvdwOfSyntyEYBPgq8LiKOB04DjgP+PqU0\nllK6nWoIONQFKaVHUkqjqerKlNJOgJTSDcC1VH/h1vpcSmkopbQBuBu4MaW0IaU0ClzJISMRE7g8\n+xz2AKdSDRj7nQPckFL6flbDL4FvAWdl7WcBX0kpbUopDQGfBKY7pbMxpfSNlNK+lNK/A9dRDVWk\nlG5KKd2dvd4MrJigv19LKW1NKY1QDW+Ls+Wvohoqzs8+41uZ+DOWWpLhQiqvT6aUjqv583RqTovU\nSildD/wt1X8Zd0XE9yLihCn2/TxgS832FWA4Wz4fqKSU9tWsv3WCfWyrfRMRZ0fEmuyUxB7gD6n+\ni/3AYYDOmvcDwCOHvJ9uJOGs7HN4CTAHqJ3QeSLw5kPC2F8Cx2ftz+HgPndS7fNUtkzw/rkAEfHK\niPhpdpplL/C/OLi/HNK/fuDY7PV8YGdKaXyKY0kty3AhtYmU0kUppdOABcAI8E9Z074JVt9O9Zcx\nANmIxROy5TuB4yOi9u+HBRPs48B+I+L5VEcJPg48MwsA1zBJGDoM++dcbAI+AVxUM89hG7CqNoyl\nlOamlP44a98JHAhc2fyOJzK1QwPaAmBH9voK4OfAidmpkc8w8/7uBJ5zyGTYiT5jqSUZLqQ2EBGv\niIjTIuLxVP813g+MZc2PACcc8ovsMuDTEfG8iDgGuBD4aUrpYapzDPYCn8omfr4SeMc0JewfcejK\n6nkT8PpDy6yzexNKKV0FdAMfzhZdCpwREWdmdT8+Il4eEa/I2lcCH4qIkyPiSVTnoow/ds8HOTki\n/ltEHBURbwZ+H/jXrO1YYG9KaSgiTgE+MIvybwV2A/8zq/O3qc7NkNqC4UIqp8nmVky2/KnA16n+\nst1Jddj9I1nb94BHqZ4u2Z0t+weqcyJuAR6kOvfgHICU0hjwJ8AfU/0F+Hmqv7hrTyEcVEdKaT3V\n0wLXZTW8HfjBNLXP9rLaidb/PPCJiHhSNt/jDcC5QCX781V+dSriYqqh6hfA/UAH0DvNMX8M/A7V\nz2EZcHZK6cGs7dzs2I9S/exXzqDeakP1M34L1Qmju6gGnW9OU4vUMiKl0l02L6lkIuIiqn9fnFt0\nLZLKz5ELSY8REa/OTplERLyW6pUW3y26LkmtoenhIiK+nN1c5sBtjLPlL4qIm7Ib1dwWES+ts+2U\nZvdJakMvpDovoBf4GvA/Uko/K7YkSa2iiJGL71G94czmQ5YvBy5KKb0E+ALVmef1tF3SkKqlI0hK\n6dsppeellI5JKf1aSumfi65JUusobM5FRDwE/ElKaW1EzAM2Acftv7Y+IipUQ0hvPW01k64kSVIT\nlWXOxfOZ+KY9Cw6jTZIkFWBO0QU0W0Q8g+rlapuBoWKrkSSppRxN9QZ816aUdk22UlnCxTZgfkQ8\nrmYUYgHVUYjeOtsm8wayBzRJkqS6nE31LrUTKkW4SCl1RUQH1Zv4XBIRfwps2z9vot62SWwGuOyy\nyzjllPa6sGTp0qUsW7as6DIaol371q79gvbtW7v2C9q3b+3aL2h+39avX8+73/1ueOxFGQdperjI\nbsbzZuDZwLUR0ZtSOhn4C+BbEfFpoIfq45H3q7dtIkMAp5xyCosWTfcAxtYyd+7ctuvTfu3at3bt\nF7Rv39q1X9C+fWvXfkGhfZtyWkHTw0VK6S8mWX4f1ccQ59YmSZKaryxXi0iSpDZhuJAkSbkyXJRE\npbfCedefR6W3Uvc+lixZkmNF5dKufWvXfkH79q1d+wXt27d27ReUt29H3FNRI2IRsGbNmjWlmuDT\nUelg8YrFrHn/GhbNL09dkiTt19HRweLFiwEWp5Q6JlvPkQtJkpQrw4UkScqV4UKSJOXKcCFJknJl\nuJAkSbkyXJRApbfC8tXLiy5DkqRcGC5KoNJXYUXHCgCWr15+WPe6kCSpaIaLklnRsYJP/PQTBgxJ\nUssyXJTQ5esup9JnuJAktSbDhSRJypXhoqS6BrqKLkGSpLoYLkqqe6C76BIkSaqL4aIEuvodpZAk\ntQ/DRQls2LWh6BIkScqN4aIEeod7iy5BkqTcGC4kSVKuDBeSJClXhosSmOi0yOY9m5tfiCRJOTBc\nlMDDfQ8/ZtndD99dQCWSJB0+w0XBKr0V7th5x4TLJUlqRYaLglX6KqzvXv+Y5cPjwwVUI0nS4TNc\nlNTo+GjRJUiSVBfDRUkNjA4UXYIkSXUxXBRsslt/jyZHLiRJrclwUbDuQR9QJklqL4aLktrWs427\nKncVXYYkSbNmuCipsTTGjdtuLLoMSZJmzXBRYnuH9hZdgiRJs2a4kCRJuTJcSJKkXBkuJElSrgwX\nkiQpV4aLgu0Z3FN0CZIk5cpwUbCprgjxyaiSpFZkuCixLXu3FF2CJEmzZrgosUf6Him6BEmSZs1w\nUWKrH17NHTvuKLoMSZJmxXBRsN7h3inbb9txW5MqkSQpH4aLgvWOTB0uvAW4JKnVGC4K1j/SX3QJ\nkiTlynBRsF0Du6Zsn+60iSRJZWO4KNjQ2NCU7dOdNpEkqWwMF5IkKVeGi4INjw0XXYIkSbkyXBRs\neNxwIUlqL4aLgo2OjxZdgiRJuTJcFGxk38iU7V6qKklqNYaLkrtj5x0+HVWS1FIMFwUb3Tf1aZH1\n3eup9BkuJEmto3ThIiLeFBFrIuLOiFgbEX+WLZ8XEddExH3Z8tfUbDNpW5lVeit09XcVXYYkSbma\nU3QBE7gU+N2U0r0RcQKwISK+D1wA3JJS+sOIeAVwZUScmFIaBz4/RVtpre1cy95hnx0iSWovpRu5\nAPYBT89ezwW6gRHg7cBFACml1cBO4Pey9d5xSNuOmrbS2rRrU9ElSJKUuzKOXLyL6shDP/A04Ezg\nWGBOSqmzZr3NwIKIOG6Cti3AgibVWzefeCpJakelGrmIiKOAvwHemlI6EfgD4DKqISgKLK0hfCiZ\nJKkdlW3k4lRgfkrpJqie4oiI7cBCYDQinlUzQnEisCWltDsixiZo2zrVgZYuXcrcuXMPWrZkyRKW\nLFmSX2+msbN354zW6xpw0qckqblWrlzJypUrD1rW09Mzo23LFi62AfMj4tdSShsi4kXAC4ENwPeA\nDwDnR8QrgecA/5ltN1HbDVMdaNmyZSxatKhB3ZiZPYN7ZrRe90B3gyuRJOlgE/2Du6Ojg8WLF0+7\nbanCRUqpMyLeD3w3Isapnrb5UEppe0R8Erg0Iu4DhoGza64GmaqttKZ73LokSa2oVOECIKX0HeA7\nEyzvBN4wyTaTtpWZT0SVJLWjUk3oPNL4RFRJUjsyXEiSpFwZLiRJUq4MFwUaGZ/6ceuSJLUiw0WB\npnsi6n4zvWRVkqQyMFy0AG8TLklqJYaLFuBtwiVJrcRwUaCx8bEZrdc7YriQJLUOw0WB+kf7iy5B\nkqTcGS4KNDg2OKP1uvt9togkqXUYLgq0L+2b0Xq7B3c3uBJJkvJjuGgBPoNEktRKDBcFmunIhc8g\nkSS1EsNFQSq9FQZGB2a07uj4zG62JUlSGRguClLpqzCWZnYpqleVSJJaieGiIF39XTNe1zkXkqRW\nYrgoSPfgzC8v3TW4i0pvpYHVSJKUH8NFC+gf66fSZ7iQJLUGw0VB1j28rugSJElqCMNFQdZ2ri26\nBEmSGsJwUZDeIR9GJklqT4aLgnhjLElSuzJcSJKkXBkuJElSrgwXkiQpV4aLFjGbO3pKklQkw0WL\nmM0dPSVJKpLhogCV3go7enfMapvuAcOFJKk1GC4KUOmrsLNv56y22fno7NaXJKkohosW0TviTbck\nSa3BcNEi+kf6iy5BkqQZMVwUoJ4rP3YN7GpAJZIk5c9wUYB6rvwYGhtqQCWSJOXPcCFJknJluCiA\nl5VKktqZ4aIAvcNe+SFJal+GC0mSlCvDhSRJypXhokUMjw0XXYIkSTNiuGgRPhVVktQqDBctYteQ\nN9GSJLUGw0UB6rlapHuwm609WxtQjSRJ+TJcFODhvofr2s77Y0iSWoHhogADowNFlyBJUsMYLgow\nODpYdAmSJDWM4aIA23u217WdV4xIklqB4aLJtvZs5e6uu+vatp6nqUqS1GyGiyZzUqYkqd0ZLlrI\nnsE9RZcgSdK0DBctZO/Q3qJLkCRpWoYLSZKUK8NFC6nnzp6SJDWb4aKF1HtnT0mSmql04SIinhAR\nX4mI+yLi7oj4drb8RRFxU0RsjIjbIuKlNdsc2nZKcT1oHO/sKUlqBaULF8AFwL6U0skppZcDH8+W\nLwcuSim9BPgC8K2abQ5tu6SJ9TaNd/aUJLWCUoWLiHgy8F7gM/uXpZQ6I2IesBi4PFv2feD5EfHC\nqdqaXX+j9Qz1FF2CJEnTKlW4AE4CdgOfiYg7IuKGiDgDeD5QSSntq1l3K7BgmrbSOZxbeO8e2p1j\nJZIkNUbZwsUc4ATgnpTSK4GPAN/JlkeRheXlcG7h7dUikqRWMKfoAg6xFRgHrgBIKd0VEQ9RDRzH\nR8TjakYoFmTr9wLzJ2mb1NKlS5k7d+5By5YsWcKSJUty60zextN40SVIko4QK1euZOXKlQct6+mZ\n2en5UoWLlNKuiPgZ8Ebgmoh4AXAicCPQAZwDXBIRfwpsSyk9CBARk7ZNZtmyZSxatKhxnWmAzv5O\n7qrcxanzTy26FElSm5voH9wdHR0sXrx42m3LdloE4APAJyJiLbAKeH9KqQL8BXBuRGwE/hp4T802\nU7WVypY9W+rediyNcW/3vTlWI0lS/ko1cgGQUnoIOGOC5fcBr5pkm0nbymbjro1FlyBJUkOVceSi\nrXkjLElSuzNctBgfuy5JKjvDRYvxseuSpLIzXDSZt/CWJLU7w0WTDY0NHdb23khLklR2hosmGx4b\nPqztfey6JKnsDBdNNjx+eOFi18CunCqRJKkxDBdNNjo+eljbH+5pFUmSGs1w0WR7Rw7vao/DPa0i\nSVKjGS6aaGvPVjbv3XxY+9g96GPXJUnlZrhoou6B+h+3vl/nQGcOlUiS1DiGixYzNOqcC0lSuRku\nWsw440WXIEnSlGYdLiLiqpkskyRJR6Z6Ri4WTLDshYdbiCRJag9zZrpiRJwL/AVwckR01DTNBe7N\nuzBNLKVUdAmSJE1pxuEC+DGwEfhnYGnN8keBtXkWpcntS/uKLkGSpCnNOFyklLYAW4BTGleOppNw\n5EKSVG6zGbkAICJOBP4HcFLt9imlM3KrSpIktaxZhwvgu8DPgK+C10U2m3MuJEllV0+4ODql9Knc\nK9GMjKdxKr0V5h87v+hSJEmaUD2Xot4TERNdjqpp3LLtlsPeRyJR6avkUI0kSY1Rz8jFPODuiLgF\nOHAv6pTSmblV1aZu33F70SVIktRw9YSLy7I/mqWB0YGiS5AkqeFmHS5SSpc0ohDNXFd/V9ElSJI0\nqXouRb14ouUppfcefjmaie7Bw390uyRJjVLPaZE1Na+PBt4GdEyyrhpgz+CeokuQJGlS9ZwW+Vrt\n+4j4Z+B7uVXUxvYO7s1lP+u71ueyH0mSGqGeS1EPNQS8IIf9tL2e4Z5c9tPZ35nLfiRJaoR65lxc\nWPP2KOAVwC9zq6iNjY6P5rKfwdHBXPYjSVIj1DPnovaf32PAPwGr8imnvY3sG8llP4/0P5LLfiRJ\naoR65lyc34hC2l2lt5Lb6YzOPk+LSJLKa9ZzLiLi2Ij4WkTcl/35akQc24ji2kmlr5LbJaQDY96M\nS5JUXvVM6Pw61RGPdwBvz/bx9TyL0tR8MqokqczqmXOxMKX08pr3H4yIu/MqSNMbTz7pXpJUXvWM\nXBxVexoke31UfiW1pzxv2b0v7cttX5Ik5a2ekYtLgFsj4jvZ+3cA/5JfSe0pz1t2Gy4kSWU243AR\nEU8FjkspfTEi7gFemzV9HZ+SOq08b9ntnAtJUpnN5rTIF4DFACmla1JKH08pfRzoBS5oRHHtZPPe\nzbnta3h8mEpvJbf9SZKUp9mEi99NKX3/0IUppUuB382vpPaU55yLsTTG2kfW5rY/SZLyNJtwMdUl\nCk4CmMbAaL73pli9c3Wu+5MkKS+zCRdzsnkXB4mIucDj8ytJM7Gjd0fRJUiSNKHZhIt/BS6NiKfv\nX5C9/pesTU3UP9JfdAmSJE1oNuHi74C9wLaIuDMi7gS2UZ3Q+blGFNdO8n6S6Y5HHbmQJJXTjC9F\nTSmNA38eEZ8FFmWLO1JKDzSksjbTM9Qz/Uqz8Ojwo7nuT5KkvNTzVNQHAAPFLA2PDxddgiRJTVHP\n7b9Vh9Hx0Vz3NzI+kuv+JEnKi+GiSUb25RsGtj26zRtpSZJKyXDRJGPjY7nub/fQbip9hgtJUvkY\nLppkNOV7WkSSpLIyXDTJyJhzJCRJRwbDRZN4tYgk6UhhuGiCSm+FvpG+3Peb58PQJEnKSynDRUS8\nJyL2RcRbsvfzIuKaiLgvItZGxGtq1p20rSwqfRWGxody32/3YHfu+5Qk6XCVLlxExAnA+4BbahZ/\nHrglpXQy8F7giog4agZtbW3P4J6iS5Ak6TFKFS4iIoBvAB8GamdAvgO4CCCltBrYCfzeJG07atra\n2vqu9UWXIEnSY5QqXAAfBX6RUrpz/4KIOA6Yk1LqrFlvM7BgkrYtwIJmFFu0LXu3FF2CJEmPMetn\nizRKRLwMeBtQujkTZdXZ3zn9SpIkNVlpwgXVUHECsCk7PXI8sAI4DxiLiGfVjFCcCGxJKe2OiIna\ntk53sKVLlzJ37tyDli1ZsoQlS5bk0ZeDdA005qqO/tH+huxXkqSVK1eycuXKg5b19MzsCd+lCRcp\npYvI5k4ARMR1wIUppX+LiN8CPgCcHxGvBJ4D/Ge26vcmaLthuuMtW7aMRYsWTbdaLroHvKpDktRa\nJvoHd0dHB4sXL55229KEiwkkILLXnwQujYj7gGHg7JTS+Aza2trQWP6Xt0qSdLhKGy5SSmfUvO4E\n3jDJepO2lUWjLhnt7O+k0lth/rHzG7J/SZLqUbarRdrSvV33NmS/faN9PhlVklQ6hosmuL/7/obt\nu1GTRSVJqpfhogkaObrgZFFJUtkYLppgYGyg6BIkSWoaw0UTjI837uKVB3Y/0LB9S5JUD8NFE4yM\nj0y/Up0e7nu4YfuWJKkehosmGEtjRZcgSVLTGC4arNJboW+kr+gyJElqGsNFg63tXMvQeOPupNnd\n79UikqRyMVw02IbuDQ3df6XXm2hJksrFcNFgvcO9Dd3/8PhwQ/cvSdJsGS4arNHhopFXokiSVA/D\nRYP1jjQ2XIzuG23o/iVJmi3DRYvrH+4vugRJkg5iuGiw/pHG/vLfM9yYx7lLklQvw0WDDYw29rki\n/SP9XjEiSSoVw0WLG2ectZ1riy5DkqQDDBcNtndwb8OP8dnrP+vohSSpNAwXDdYz3NPwY9y8/WYq\nfYYLSVI5GC4arNFzLiRJKhvDRYMNjDUnXNy67damHEeSpOkYLtrEPV33FF2CJEmA4aLhxsfHiy5B\nkqSmMlw0WLOe/dHom3VJkjRThosGqvRW6Blp/NUi4MRRSVJ5GC4aqNJXYXBssCnHGhxtznEkSZqO\n4aKBuvq7mnasnY/ubNqxJEmaiuGigTbt3tS0Yz3U85B36ZQklYLhooH2DjX+1t8HjjW817t0SpJK\nwXDRQNsf3V50CZIkNZ3hooG29xguJElHHsNFAw2NDRVdgiRJTWe4aCMbujcUXYIkSYaLRhoeG27q\n8R7c82BTjydJ0kQMFw3UN9LX1ON5KaokqQwMFw3UrFt/7/eTB35iwJAkFc5w0UDNft7H/Xvu914X\nkqTCGS4aaHyfj1uXJB15DBcNUumt0DPU3NMikiSVgeGiQSp9FUbTaNOP6+WokqSiGS7azLKblzmp\nU5JUKMNFgxQ1sXL1w6ud1ClJKpThokGa+URUSZLKxHDRhroGuoouQZJ0BDNcNMiewT2FHXvTrk2F\nHVuSJMNFgxR5WsRTMpKkIhku2tA1m67xihFJUmEMFw2y9uG1hR375u03e8WIJKkwhosGqPRWuHHb\njUWXIUlSIQwXDVDpqxQ+cuAVI5Kkohgu2tSK1SucdyFJKoThok2t2rCq8NETSdKRyXAhSZJyVapw\nERFPjIgrI2JDRNwZEddGxElZ27yIuCYi7ouItRHxmprtJm0rgvMdJElHslKFi8zylNKvpZR+E/gh\n8I1s+QXALSmlk4H3AldExFFZ2+enaGu67oHuog4tSVLhShUuUkrDKaUf1yy6FTghe/124KJsvdXA\nTuD3srZ3HNK2o6at6Yq89Xet5auXO6lTktR0pQoXE/gIcFVEHAfMSSl11rRtBhZM0rYFWNC8Mg9W\nlttvr+hY4aROSVLTzSm6gMlExKeBk4D3A08uuBxJkjRDpQwXEfFx4K3Aa1NKQ8BQRIxFxLNqRihO\nBLaklHZP0rZ1qmMsXbqUuXPnHrRsyZIlLFmy5LDrL9OpCCeXSpLqsXLlSlauXHnQsp6enhltGyml\nRtRUt4j4KHAW1WDRU7P8Yqph4vyIeCWwCjgxpTQ+VdsE+18ErFmzZg2LFi1qSB/+6PI/4ur7r27I\nvmfrK3/4FboHujl38bnMP3Z+0eVIklpYR0cHixcvBlicUuqYbL1SzbmIiOcCXwLmAtdll6PekjV/\nEnhVRNwHXAycXRMepmprqkpvhdU7Vxdx6AntHdrL+Tec79wLSVLTlOq0SEppB5MEnuyUxxtm29Zs\nlb4Kjww8UnQZB/xg/Q+KLkGSdIQp1chFO+jqL9cch9UPV0dRNnZv5LzrzyvVfBBJUnsyXORs0+5N\nRZcwoQf2PODpEUlSUxguclaWe1wcyhELSVKzGC5ytv3R7UWXMKEte7cUXYIk6QhhuMjR1p6t/J81\n/6foMia049EdQPnmhEiS2o/hIkfdA93sY1/RZUzors67gPLOCZEktQ/DxRGmrHNCJEntw3BxhKn0\nVrwkVZLUUIaLHLXCfIZ1D6/zklRJUkMZLnLUPdhddAnT2jW0q+gSJEltznCRo817NhddwrQe3PNg\n0SVIktqc4SJHdz18V9ElTGtofAiAZbcsc96FJKkhDBc5qfRWuHX7rUWXMWOXrbvMeReSpIYwXOSk\n0ldhe285784pSVIzGS4kSVKuDBdHsAtvudB5F5Kk3BkujmCXr7ucjoc7ii5DktRmDBdHuDO/cyZb\ne7YWXYYkqY0YLnJy7aZriy6hLiPjI9y09aaiy5AktRHDRU5u2ta6v6Af2PNA0SVIktqI4SIHld4K\nax9ZW3QZdbt/1/2c/q3TuatS/puASZLKz3CRg0pfhW2924ouo243bL2BG7bcwL3d9xZdiiSpDRgu\nDlOlt8Ly1cuLLuOwbN67GYA9g3uKLUSS1BYMF4ep0ldhRceKosvIxfqu9UWXIElqA4YLHfCTB37i\nTbUkSYfNcKED7t9zvw8zkyQdNsPFYerq7yq6hFxdeMuFfOzajzmCIUmq25yiC2hlld4KK9a0x3yL\n/S5fdzkAZy88m/nHzi+4GklSK3Lk4jBU+iqs2rCq6DIaYvnq5Y5eSJLqYrjQhFZ0rHD+hSSpLoYL\nSZKUK8NFnSq9Ff7m539TdBkNteyWZZ4akSTNmuGiTpW+Ctfcf03RZTTUZesu89SIJGnWDBeSJClX\nhos63brt1qJLkCSplAwXdaj0Vli+prUfVjZTZ1xyBgu/vtDHsUuSZsxwUYdKX4W1nWuLLqMpeoZ7\nWNe1jqs3XV10KZKkFmG40Ix8997veuWIJGlGDBd16Bpor+eJzMTazrV8+EcfNmBIkqZluJilSm+F\nz9/4+aLLKMSqDav44NUfNGBIkqZkuJiltZ1ruX7z9UWXUZirNl7Fku8vMWBIkiZluJil7oHuokso\n3A1bbuC0b57GTx/4Keddf55BQ5J0EB+5Pkt3V+4uuoRS2NKzhbP+71l0D3Xzlpe8xcezS5IOcORi\nFr56+1f54i1fLLqM0ugeqo7i/N0Nf8ddlbv42LUf42PXfsyRDEk6wjlyMQs/3vTjoksopSs3Xsm6\nznXcv+co9iAaAAAMlUlEQVR+AF5/0usdyZCkI5gjF7PQ1X/kXYI6U/uDBcDZq87mpw/8tMBqJElF\nMlzMUKW3wj2d9xRdRkvYNbiL23bcVnQZkqSCeFpkhv7lzn9hYHyg6DJaxldu+wobuzfyrKc8i4+/\n6uMzPk2ytWcrAAvmLmhkeZKkBjJczMA1m67hc//5uaLLaCmdA51ctu4yAF7/otcDsHzNcs5dfO5B\nQWP/Q+DOXXwuo/tGeclXXwLAxg9vNGBIUosyXEzjmk3X8KYr3lR0GS3tA//+AfYO7mXP8B5OfsbJ\nnPUbZx1oW9u5lvNvOJ/Tnn8a8548j6GxIaB6P5HJwoWjG5JUbs65mEKlt8In/+OTRZfR8h7a+xB7\nhvcA1cme874wj1O+egqLly9m9Y7VAGzatYkLb7nwwDbLVy+f8JLWO3bewQu//EJe/E8vZumPl3rZ\nqySVkCMXk6j0Vvjwjz58xDxavZm6B7vpHqzeI2P7o9sB+N83/28292w+sM6KjhWc+dIzAfjSzV+i\nb6QPgG2PbmM8jTOexvnH2/6RgdEBzjv9vElPtXhJrCQ1n+FiEms717Jqw6qiy5iddcBvFF3E7HQO\ndAIcFCz2e98P3scxjz+GDbs3TNq3FR0r6B7o5plPfuaBZf2j/Vy+7nJOe95pfOnmLwFwzsJzuGrj\nVRMGjkpv5cB6s5l8moeVK1eyZMmSph2vmdq1b+3aL2jfvrVrv6C8fWubcBERLwIuAZ4J7AX+a0pp\nfT37qvRWOPv/np1nec3RguFiKtt7t//qzRR9mywEXnDTBVy3+ToAbtt+Gzdtv4nTnncaj/Q9wgeu\n/gALn72Q804/j0pfhQtvrZ6SOXvh2YzuGwWaM6ejrH8x5KFd+9au/YL27Vu79gvK27e2CRfAcuCi\nlNKlEfE2qkHjt2azg0pvhfOuP48fbvwhu4Z2NaRINc/+YAFw0/abAHj7997OMY8/hkp/hVt33Er3\nQDdrH/nVqa/PXv9ZfnT/j0gkVr1zFccfczxX33e1p1gkaRbaIlxExDxgMfA6gJTS9yPiqxHxwpTS\ngzPdz3889B+s6FjRqDJVAr0jvfSO9B54f+ioxw/u+8GB129Z+RaCIJHYtGsT7/r1d/Hy418+4YhG\n7TwPmPiyW0k6UrRFuACeD1RSSvtqlm0FFgCThoutPVt5pP8RLvjFBVy14SrGGW90nWoxiQTAFfdc\nwRX3XAHASXNP4q9e9VdcfOfFvOBpL2Dz3s2c+LQTWbVhFbdvv50H9z7Ixl0bqyNhNZNNJ5poOjg6\neODSWjj4VEyZLrktUy2SijPTx2C0S7iYjaMBfucffofR40aLriVfQ8DOootokBL17YGdD/CX6/8S\ngDu5E4AOOgC4Zuc1B9ZbsXMFK/59BU846gmMjI/wOB7HPvZxwfcuYF/axxOOegJ9v+zjhL8+4cA2\n84+ZT/9oP0fPOZrO/upk12OfcCwLn72QLXu30DXYxbwnzePFz3wxld4KL533Ul77gtfyldu/wti+\nMRY+eyHvX/z+6vHXrOCXXb/kOcc+h529O5l/7Hy6+7v51Ks/xXFPPo5L774UgDe/+M1cv+V6Tj/h\ndK7edDWDY4MH9TelxLUPXMvj4nF8/Y++zo1bbuTMU85k3lPm0dXfxaV3X3rQNk+a8yTOefk5dO7q\n5NwV53L6Cadz/ZbrOfOU6tU/tce9etPVB70eHBvkSXOedFDbOS8/h3lPmTfhd7H/+PvXA1i1ftWB\n+vavc+iyqUy3fk9PDx0dHdPu53DNpO6J1jmc/tb27dDPdib7KqtGfGez/ZwbpVk/j/vdcfcd+18e\nPdV6kVJqfDUNlp0W2QQct3/0IiIqwH859LRIRJwFXN78KiVJahtnp5SumKyxLUYuUkpdEdEBnANc\nEhF/CmybZL7FtcDZwGaq/x6WJEkzczRwItXfpZNqi5ELgIg4GfgW8AygB3hPSuneQouSJOkI1Dbh\nQpIklYPPFpEkSbk6osJFRLwoIm6KiI0RcVtEnFJ0TZOJiC9HxEMRsS8iFtYsP7QPL62zrZC+R8QT\nI+LKiNgQEXdGxLURcVLWNi8iromI+yJibUS8pma7utqaLevPXVnfboiIU7PlLf291dTznuxn8i3Z\n+3b4zjZHxPrsO+uIiLdny1v6O4uIJ0TEV7LP+O6I+Hab9Ou4mu+qI6tnJCKe1uo/jxHxpohYk/Vv\nbUT82eHUXmi/UkpHzB/gZ8A52eu3AbcXXdMUtb4aeA7V+3QsnEkf6m1rcr+eCLyx5v2HgOuy1xcD\n/3/2+hXANuCo7P0362kroH9PrXn9VuCudvjesuOfANyU/XnL4XwvJfvOHgR+Y4LlLf2dAcuAL9e8\nf1Y79GuCfn4M+EH2uqX/DgF2AS/LXp8ADAJPacV+Ff6D0cQvbR7VZ448rmZZBXhh0bVNU/dDZOFi\nqj7U21aC/i0GHsxe9+7/CzB7fxtwRh1tt+5vK7hv/xXoyD7/nlb+3oAAfgr8JnAdvwoXLf+d1f4/\nVrOspf9fA56c/cwd0079mqSvvwT+uB1+HoEu4NXZ64VUw8DjW7FfbXEp6gzVdRfPkpmqD4/W2VZ0\n3z8CXBURxwFzUkqdNW2bgQV1tG2h2rdCRMQlwO8DCXgT1e9tZ4t/bx8FfpFSujMigOrwNG3ynQHf\njmrHbgc+Rev/v3YSsBv4TET8ATAAnE81ILRyvw4SEa8CngZc3SY/j+8CroyIfqr9OhM4doL6NlPy\nfh1Rcy5ULhHxaap/CX666FrylFL685TSAuBvgC9ki6PAkg5LRLyM6jD4/yq6lgZ5TUrpVGAR1WHp\nS7LlLfudUb2H0QnAPSmlV1IN8d/Jlrdyvw71XuDbh4SelhQRR1H9O+OtKaUTgT8ALqNFv7MjKVxs\nA+ZHRG2fF1BN361iqj7U21aIiPg41TkJb0wpDaWUdgNjEfGsmtVOBLbU2Vb495pSupTqCEarf2+v\nofqLalNEPAT8DrACeAdt8J2llLZn/x0H/pHqfKettPZ3thUYB64ASCndRfX0zwnA8S3crwMi4ilU\nfwYvBmiDv0NOBeanlG4CSCmtBrZTPT0y2mr9OmLCRUqpi+r573MAYuq7eJbSVH2ot62AbhARH6U6\n/Pe6lFJvTdP3gA9k67yS6oTW/6yz7YbG9uKxImJuRMyvef9WoDv7/NfQot9bSumilNJzU0ovTCm9\ngOp52/8vpXQRrf+dPTki5tYsOgvoSCl109rf2S6qEzDfmNXxAqq/WG6sp/ay9OsQ76I6Yfq+mmWt\n/PO4P8D9WlbHi6jOd9kwSX3l7lczJnaU5Q9wMnAzsJHqudWXFV3TFLVeRPWHbYTqxKn7putDvW1N\n7tdzgX1UnwXTAdwJ3JK1PYvqLWXvA9YBv1uzXV1tTe7bAqqTqe4G7gJ+wq8m47b093ZIP3/OryZ0\ntvp39oLs5/Cu7Hu7EljQDt9Z1refA2uz/8/e2g79qqnnRuDPDlnW6j+P76z5vu4G3tmq/fIOnZIk\nKVdHzGkRSZLUHIYLSZKUK8OFJEnKleFCkiTlynAhSZJyZbiQJEm5MlxIkqRcGS4kSVKuDBeSDltE\nXB0RH5xg+V3ZbdAn2+5vI+LCxlYnqdkMF5Ly8E2qT6g8ICJeARwP/FshFUkqjOFCUh5+CDwvIn69\nZtl7gG8Dp0TELyJidUTcExGfnmgHEfHnEXFlzfs3R8R1Ne/fHRG3Zvu5PiIWZst/O1vWERFrI+Lc\nxnRR0kzNKboASa0vpTQWEZcC/w1YGhFPBJYAp1F9bPQZKaXRiDgauDki/iOldPtEu5rofUS8Ktvf\na7L9vJrq48R/HfgU8MWU0neydeciqVCGC0l5uRi4PiI+AbwN+GVKaWNEzAP+OSJOpfpE3OcBp1J9\nsuZM/QmwELgtIiJb9rQsxFwH/M+IOBn4eUrpppz6I6lOnhaRlIuU0nrgfuAtVE+JfCNr+nugC3h5\nSulU4Abg6Al2MQYcVfO+dp0ALkkpLUop/Wb253kppeGU0peBNwM7gb+PiK/l2jFJs2a4kJSni4FP\nA68EvpstezqwPaWUIuIlwOsm2fZ+YGFEPDEi5gBn1bT9EHh3RDwfIKoWZ69PTiltSSl9k2qQ+e3c\neyVpVjwtIilP3wGWAf+aUhrIlv0dcGlEvAvYDPxsog1TSrdFxI+AdcBDQAdZUEgp3RgRfw1cGRFH\nAU8ArgbWAB+OiDOAEaqjHx9tUN8kzVCkdOj8KUmSpPp5WkSSJOXKcCFJknJluJAkSbkyXEiSpFwZ\nLiRJUq4MF5IkKVeGC0mSlCvDhSRJypXhQpIk5cpwIUmScmW4kCRJufp/sfhraeIUuvQAAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1cbc7365320>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhcAAAF7CAYAAABhB6n0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3X+cXHV97/HXR6JyqxilwjVVQlCLoi1iIrS2ol6q1apX\nrb/aiLSV2xaqtt6gtyr1XqH+/gVqVRCLgiARrUGtSrlqo1coqGSBqAUSgfyALGRDkk1IyI9NPveP\nc5ZMlp3N7uyZOTM7r+fjMY+dOd/z4zMn85h553u+55zITCRJkqrykLoLkCRJM4vhQpIkVcpwIUmS\nKmW4kCRJlTJcSJKkShkuJElSpQwXkiSpUoYLSZJUKcOFJEmqlOFC6hER8a6I+HLddXSLiPjziLih\nTet+XkRsase6m2zvGRGxt1Pbk9rNcCF1gYhYGhF/N870vRFxLEBmfjAzT57Eut4TEUvaUWcnRcSq\niNgeEVsiYnNE/DQinj9mtnbev6DT90bwXgyaMQwXUnfruh+ciDioQ5tK4E8y81GZ+WjgQuCbEfGw\nDm1fUosMF1J3iweeFD0SVzS8/nBEDEbEcETcEhEviYhXAGcCL4uIrRGxpZx3VkR8MCJWR8Q9EbE4\nIh7bsK6nRcS15bp+UK57aUP73oh4c0T8HLgvIh4REYsiYkXZs7AyIt7cMP+R5TJvjIjbynk+FBGP\ni4j/W25naUQcPtn3D3wJOASYO+6MDb085eu3jnkPh0XEpRGxLiLujIhzI+KhE207It5Szr8uIs5q\naDiifB/rI+LeiPh2RBzZ0P7FiLig3M9bIuLmiHheQ/vsiPhqRGyKiP8EnnuA/SD1FMOF1FsSICJe\nCCwEjsvM2cALgBWZ+U3gA8C3M/OQzHxUudyZwEuA3wOOKqd9uVzXLOBbwHeAXwfeBZzKg3tNFpbb\neRSwHVgFPL/cxl8CH42IZ49Z5vnA04ETgP8JXA78HfBYYHdZ1wGVvSWnAneW2x3PeL08jdO+Bayj\neP+/DRwLvHuCzR4CPLOc/78Bp0bEKWXbQ4CPA48HjgS2AZ8fs/zrgM8Cs4FLgYsa2v6JYj/OBU4C\n/myCOqSeY7iQuseHImJjw2MTzQ+L7AYeDvx2RMzKzDsz81cTrPsNwHsz867M3A6cAbwwIh4HPBs4\nFPhAZo5k5k8pQsBYH87MezJzdxauyMx1AJn5I+AqijDR6L2ZuSMzbwFuAq7OzFsyczdwBTD/APvk\nyxGxkeLH+6PAOzNzpMm80WQ6EXE88OTM/PvM3JmZm4APAhONYQlgdP5bgU8Dp5Tvd3VmXlXui/vK\ndT1nzPLfzcwfZ2YCXwTmRsRjIuIhFMHjHzJza2beXb43acYwXEjd452ZeWjD4zE0+cHMzB8C7wHe\nCwxFxNcau+XH8QRgdcPyg8DOcvocYDAzG89WWDPOOtY2voiIkyNiWXlYYBPwRxQ9Eg9sBljf8Ho7\ncM+Y14+coGaA15f74mCKEPSxiPjDAywzniOBxzSGN+BfgMMmWGZHZt7b8Ho1RU8FEfHYiPhyRKyJ\niM3Aj4CHR8QhDfPf3fB8G8W/5SEU++ih7L+PVyPNIIYLqUdl5vmZ+WyKrvVdwKfKpvFOabwTmDf6\nouyxeFg5fR3wuPJ/1KPGG9fwwHoj4giKbv63A48tg9CVTNB70KIH1peZNwHXAC9tMu824NcaXs9p\neL4WuGdMeHt0eUipmYMbx6VQBJS7yucfAv4LxWGpR7NvzMRk3v8GYKRcX+O6pRnDcCH1oIh4VkQ8\nuxyQuJPih3X0cME9wJFjzuq4FDgzIp4QEY8EzgG+V3bJXwdsBt5VDvw8nqLbfiKjPQ5DZT0vAcb2\nKFQaNCLit4ETgeVNZhkATomIgyLiOMpDGKWfAWsj4r3l+x8ddPriCTaZwAcj4uCIeArwJor9CEUP\nxHZgS0T8OnDWZN9H2UP0VeAfy4Gdv0ER0qQZw3AhdYdmYyuaTX8UxWDBDRQ9D3OAt5ZtXwO2UBwu\n2VhO+yDFmIhrgduBg9g3fmAEeAXw34GNFP8rv4QitIxbR2beDLwfWFrW8FrgmweovZXTakfPttgC\nfAP4TGZe2GTev6UYsDo6nuKihnr3Ai+jOKxxc3ko41+BJ02w7S3AjRT764fARZn5pbLtPcBvltv6\nMcVg2ANpfP9/SxEIVwPfpzgTRpoxohhrJEn7RMT5FN8Pp9Vdi6TeY8+FJCLiOeUhk4iIPwBeT9F1\nL0lT1vFwERGfjIg7xrngzZMj4pqIuDUifhIRT2ux7ZhOvydpBngixdiLrcBngHdk5g/qLUlSr6qj\n5+JrwO/z4AvhfA44PzOfAnyE/S84M5W2i9tStTSDZeaXMvMJmfnIzHxqZp5Xd02SeldtYy4i4g7g\nFZm5PCIOA1YCh46eax8RgxQhZGsrbZl5e6ffkyRJ6p4xF0cw/kV85k6jTZIk1WBW3QV0WnlO+oso\nDsvsqLcaSZJ6ysEUF+S7aswVbPfTLeFiLTAnIh7S0Asxl6IXYmuLbc28iPKGTZIkqSUnA5c1a+yK\ncJGZQxExQHFRn4sj4jXA2tFxE622NbEK4NJLL+WYYzyxpB0WLVrEueeeW3cZM5r7uP3cx+3l/m2/\nduzjm2++mTe84Q3Q/O7EQA3horw4z0uB/wpcFRFbM/No4HTgoog4ExgG3tiwWKtt49kBcMwxxzB/\n/oFuyKhWzJ49233bZu7j9nMft5f7t/3avI8nHFbQ8XCRmac3mb6C4tK9lbVJkqTO65azRbre4CCc\ndVbxV5IkNWe4mKTBQTj7bMOFJEkHYrhQ5RYuXFh3CTOe+7j93Mft5f5tvzr3cd/dFTUi5gPLli1b\nNqWBLgMDsGABLFsGjkGSJPWjgYEBFixYALAgMweazWfPxRQNDTn2QpKkiRgupuiCCxx7IUnSRAwX\nU7RkSfF3aKjeOiRJ6laGixZt2FB3BZIkdSfDhSRJqpThQpIkVcpwMUmOsZAkaXIMF5M0dozFpk31\n1CFJUrczXIxjvPuIjA0Xa9d2tCRJknqG4WKM0WAx9loWW7fuP9/tt3e0LEmSeobhYozBweJCWWON\nDRd3392ZeiRJ6jWGi0kaGy527qynDkmSup3hYhIGB+Gaa/aftn17PbVIktTtDBeTMDgIy5fvP82z\nRSRJGp/hYhLGu8aFh0UkSRqf4WISxruPyMhI5+uQJKkXGC5atGdP3RVIktSdDBdjjHcIZNWqB0/b\nu7ftpUiS1JMMF2OsXPngaXfd9eBpO3fuf5EtSZJUMFyMsXnz5Obbs+fBZ5BIkiTDxYRGD5Fs2zZ+\n+3i9HJIk9TvDxRiNV+IcDQ/NLpjlzcskSXoww8UYjeFi7driJmbNLpi1bl1HSpIkqacYLiZw993F\n3VG3bBm/3at0SpL0YIaLCdx7b/F39+7x23fs6FwtkiT1CsPFBK69tvi7a1e9dUiS1EsMFxPYuLH4\n66W+JUmaPMPFJKxfP/50b14mSdKDGS4modmFtQwXkiQ9mOFCkiRVynAhSZIqZbiQJEmVMlxIkqRK\nGS7GaHaTMkmSNDmGizGa3aRsPHfeCYOD7atFkqReZLiYhsFBw4UkSWMZLhoMDu675LckSWqN4aLB\n8uVw1111VyFJUm8zXDTYsKHuCiRJ6n2GC0mSVCnDhSRJqpThQpIkVcpwMU1DQ3VXIElSdzFcTNPK\nlXVXIElSdzFcTNPatXVXIElSdzFcTNPWrXVXIElSdzFcNNi0qe4KJEnqfYaLBqtW1V2BJEm9z3DR\nwDM/JEmavq4LFxHxkohYFhE3RMTyiPizcvphEXFlRKwop5/YsEzTtqmYyu3WJUnS+GbVXcA4LgGe\nm5m/jIgjgVsi4uvAh4FrM/OPIuJZwBURMS8z9wAfmqBNkiR1UNf1XAB7gceUz2cDG4BdwGuB8wEy\n83pgHfC8cr7XjWm7q6GtrbZt68RWJEnqHd3Yc/GnFD0P24BHA68CDgFmZeb6hvlWAXMj4tBx2lYD\ncztR7L33dmIrkiT1jq7quYiIg4B3A6/MzHnAC4BLKUJQ1FhaU8PDdVcgSVJ36baei+OAOZl5DRSH\nOCLiTuBYYHdEHN7QQzEPWJ2ZGyNiZJy2NRNtaNGiRcyePXu/aWvXLgQWTqlgr40hSZqJFi9ezOLF\ni/ebNjzJ/1FHZrajppZExOHASuB3MvOWiHgycB1F6PhHijBxdkQcDywB5mXmnoj4QrO2cbYxH1i2\nbNky5s+fv1/biSfC1VdPreajjoLbb5/yW5UkqecMDAywYMECgAWZOdBsvq7qucjM9RHx18BXI2IP\nxWGbN2fmnRHxTuCSiFgB7ARObggPE7VN2s6d1bwPSZL6WVeFC4DMvBy4fJzp64EXNVmmadtU7N49\n3TVIkqSuGtBZt1276q5AkqTeZ7iQJEmVMlxM0x6vASpJ0n4MFw1GRjqzjCRJM5nhooEDOiVJmj7D\nRQN7ISRJmj7DRQPHT0iSNH2Gi2nau7fuCiRJ6i6Gi2lynIYkSfszXEyT4UKSpP0ZLhq0MubCwyKS\nJO3PcNHAs0UkSZo+w0WDVnouPMNEkqT9GS4aZNZdgSRJvc9w0cAxF5IkTZ/hYprs7ZAkaX+Giwat\n9FwYLiRJ2p/hooGHOCRJmj7DhSRJqpThYpo8LCJJ0v4MF9O0Zw8MDtZdhSRJ3cNwMU2ZsHx53VVI\nktQ9DBcNWj3EsWFDtXVIktTLDBcNHD8hSdL0GS4qsGlT3RVIktQ9DBcV2Ly57gokSeoehosGrR4W\n8WwRSZL2MVxUYP36uiuQJKl7GC4kSVKlDBeSJKlShgtJklQpw0WDVu+Kev/91dYhSVIvM1w0aPVs\nkeHhauuQJKmXGS4qsHVr3RVIktQ9DBcNWu252Lat2jokSeplhovSdC6E9atfwZo11dUiSVIvM1yU\npnuVTe+MKklSwXAhSZIqZbiQJEmVMlyUhobqrkCSpJnBcFFaubLuCiRJmhkMF6XNm+uuQJKkmcFw\nIUmSKmW4KHmVTUmSqmG4KN199/SWd0CoJEkFw0Vp+/bpLe9FtCRJKhguKrJpU90VSJLUHQwXFfFs\nE0mSCoYLSZJUKcOFJEmqlOGidP/9dVcgSdLMYLgo7dgxveW9ToYkSQXDRUUMF5IkFbouXETEwyLi\nnyJiRUTcFBFfKqc/OSKuiYhbI+InEfG0hmXGth1T3zuQJKm/dV24AD4M7M3MozPzGcDby+mfA87P\nzKcAHwEualhmbNvFHaxXkiQ16KpwERG/BpwK/MPotMxcHxGHAQuAL5fTvg4cERFPnKitk7Vv29bJ\nrUmS1L26KlwATwI2Av8QET+LiB9FxEnAEcBgZu5tmHcNMPcAbR1z772d3JokSd2r28LFLOBI4BeZ\neTzwVuDycnrUWdiBTPdsE0mSZopZdRcwxhpgD3AZQGbeGBF3UASOx0XEQxp6KOaW828F5jRpa2rR\nokXMnj37gdc33ACwsHxIktTfFi9ezOLFi/ebNjw8PKlluypcZOa9EfED4MXAlRFxFDAPuBoYAE4B\nLo6I1wBrM/N2gIho2tbMueeey/z58x94fcIJsHFjG96UJEk9aOHChSxcuP9/uAcGBliwYMEBl+2q\ncFH6G+DCiPgwRS/GX2fmYEScDlwUEWcCw8AbG5aZqG1Sdu+eXtFe50KSpELXhYvMvAM4aZzpK4Df\na7JM07bJ2r59Okt7V1RJkkZ124DO2ky352K6y0uSNFMYLiqyZ0/dFUiS1B0MFxXxrqqSJBUMFxUx\nXEiSVDBcVMTDIpIkFQwXkiSpUoYLSZJUKcNFabqHNTKrqUOSpF5nuCiNjExv+R07YM2EdzORJKk/\nGC4qtGFD3RVIklQ/w4UkSaqU4UKSJFXKcCFJkio15XAREd+YzLRes3dv3RVIkjQztNJzMXecaU+c\nbiF1m+7ZIpIkqTBrsjNGxGnA6cDRETHQ0DQb+GXVhUmSpN406XAB/BtwK3AesKhh+hZgeZVF9aqh\noborkCSpfpMOF5m5GlgNHNO+cnrbypXwohfVXYUkSfWaSs8FABExD3gH8KTG5TPzpMqq6lGbN9dd\ngSRJ9ZtyuAC+CvwA+DTgjcYbDA7WXYEkSfVrJVwcnJnvqrySGWD16rorkCSpfq2civqLiBjvdNSe\nVsVdTXfsmP46JEnqda30XBwG3BQR1wIP/Jxm5qsqq6oGVVznYufO6a9DkqRe10q4uLR8zCh7Khg9\nct9901+HJEm9bsrhIjMvbkchdaui52L79umvQ5KkXtfKqahfGG96Zp46/XLqU8WYi927p78OSZJ6\nXSuHRZY1PD8YeDUw0GReSZLUZ1o5LPKZxtcRcR7wtcoqqkkVd0X1sIgkSa2dijrWDuCoCtZTqyoO\ni9x///TXIUlSr2tlzMU5DS8PAp4F/GdlFfWwKno/JEnqda2MuRhueD4CfApYUk059ami56KK01kl\nSep1rYy5OLsdhcwEu3YV9xeZM6fuSiRJqs+Ux1xExCER8ZmIWFE+Ph0Rh7SjuE6qqufCm5dJkvpd\nKwM6P0vR4/E64LXlOj5bZVG9bGio7gokSapXK2Mujs3MZzS8flNE3FRVQb1uw4a6K5AkqV6t9Fwc\n1HgYpHx+UHUlSZKkXtZKz8XFwHURcXn5+nXAF6srqR5VjLmQJElTCBcR8Sjg0Mz8aET8AviDsumz\nzMC7pLZq06a6K5AkqV5TOSzyEWABQGZemZlvz8y3A1uBD7ejuE6qqudi7dpq1iNJUq+aSrh4bmZ+\nfezEzLwEeG51JfW2devqrkCSpHpNJVxMdP3Jnr7wdZXXptixo7p1SZLUi6YSLmaV4y72ExGzgYdW\nV1LneeErSZKqM5Vw8RXgkoh4zOiE8vkXyzbhgE5JkqYSLt4HbAbWRsQNEXEDsJZiQOd721FcL9qy\npe4KJEmq16RPRc3MPcCfR8Q/AvPLyQOZeVtbKuugX/6y7gokSZo5Wrkr6m1AzweKRqtXV7euXbuq\nW5ckSb2olct/awK7d9ddgSRJ9TJcVGxkpO4KJEmql+GiYtu21V2BJEn1MlxUzHAhSep3houKeVhE\nktTvDBcV2zPRRdIlSeoDhouK3X8/rFlTdxWSJNWnK8NFRLwxIvZGxMvL14dFxJURsSIilkfEiQ3z\nNm2ry4YNdVcgSVJ9ui5cRMSRwF8C1zZM/hBwbWYeDZwKXBYRB02iTZIkdVhXhYuICOCfgbcAjde6\nfB1wPkBmXg+sA57XpO2uhjZJktRhXRUugDOAH2fmDaMTIuJQYFZmrm+YbxUwt0nbamBuJ4ptZmio\nzq1LklSvrgkXEfF04NXA++uuZbquv77uCiRJqs+Ub1zWRicCRwIry8MjjwMuAM4CRiLi8IYeinnA\n6szcGBHjtR3wfI1FixYxe/ZsAFasGJ26sHxMz771SZLUmxYvXszixYv3mzY8PDypZSMz21HTtEXE\nUuCczPzXiPgCRZg4OyKOB5YA8zJzz0RtTdY7H1i2bNky5s8v7hz/jnfARz5SXe0vfSl8+9vVrU+S\npG4wMDDAggULABZk5kCz+bqp52KsBKJ8/k7gkohYAewETm4IDxO1TcrWrRVVXJpksJMkaUbq2nCR\nmSc1PF8PvKjJfE3b6nLffXVXIElSfbpmQOdM4s3LJEn9zHBB9WFg+/Zq1ydJUi8xXFB9GNi5s9r1\nSZLUSwwXkiSpUoaLNti8GQYH665CkqR6GC7aYGTEcCFJ6l+GC0mSVCnDhSRJqpThok28M6okqV8Z\nLoD7769+nStXVr9OSZJ6geGC9twLZO3a6tcpSVIvMFzQnoterVtX/TolSeoFhos2uewyWLOm7iok\nSeo8wwWwe3f169y7FzZsqH69kiR1O8MFsGtX3RVIkjRzGC4kSVKlDBfAjh11VyBJ0sxhuMBwIUlS\nlQwXbeRVOiVJ/chw0UaeLSJJ6keGC4rTRiVJUjUMF8DISN0VSJI0cxgu2mjJEhgcrLsKSZI6y3DR\nRoYLSVI/MlxIkqRKGS6AzLorkCRp5jBc0J5bro/yWheSpH5juAD27Gnfur3WhSSp3xgu8LCIJElV\nMlxguJAkqUqGCwwXkiRVyXCB4UKSpCoZLmjvvUU2bWrfuiVJ6kaGC9rbc3Hzze1btyRJ3chw0War\nV9ddgSRJnWW4aLPh4borkCSpswwXbbZypTcvkyT1F8NFm91zj+FCktRfDBeSJKlShosO8OZlkqR+\nYrjoAG9eJknqJ4aLDvB0VElSPzFc0P7Lf19xhYM6JUn9w3BB+8PF9dcbLiRJ/cNw0SG33lp3BZIk\ndYbhokNuu63uCiRJ6gzDRYds3Vp3BZIkdYbhokMMF5KkftH34cKBlpIkVctwYbiQJKlSfR8uOsWr\ndEqS+oXhokM2bqy7AkmSOqPvw8Uvf1l3BZIkzSxdFS4i4uERcUVE3BIRN0TEVRHxpLLtsIi4MiJW\nRMTyiDixYbmmbQfSqft+2HMhSeoXXRUuSp/LzKdm5jOBbwH/XE7/MHBtZh4NnApcFhEHlW0fmqCt\nK9x2m4NHJUn9oavCRWbuzMx/a5h0HXBk+fy1wPnlfNcD64DnlW2vG9N2V0NbV9iyxXAhSeoPXRUu\nxvFW4BsRcSgwKzPXN7StAuY2aVsNzJ3MBjp5cSvHd0iS+kHXhouIOBN4EnBmO7fTyXBx6qmwZk3n\ntidJUh1m1V3AeCLi7cArgT/IzB3AjogYiYjDG3oo5gGrM3Njk7YJf8YXLVrE7NmzWb68cerC8tEe\nIyPF9S7mTqpPRZKk+ixevJjFixfvN214eHhSy0ZmtqOmlkXEGcDrKYLFcMP0L1CEibMj4nhgCTAv\nM/dM1DbO+ucDy5YtW8b8+fN505vgvPM68tYAWLYM5s/v3PYkSarKwMAACxYsAFiQmQPN5uuqnouI\neDzwMeA2YGlEBLAjM58NvBO4JCJWADuBkxvCw0RtkiSpg7oqXGTmXTQZB1Ie8njRVNskSVJnde2A\nzk7p9D0/brmls9uTJKnT+j5cdPrKmd/8Zme3J0lSp/V9uNi5s7Pbu+OOzm5PkqRO6/twcd99nd3e\n5s2d3Z4kSZ3W9+Fi+/bObq+TF+2SJKkOfR8udu/u7PZ27Ojs9iRJ6rS+DxedZriQJM10hosO2+Ol\nvSRJM1zfh4uRkc5ur8uuti5JUuX6Plx0uidhZARuvLGz25QkqZP6PlzU4YwzYHCw7iokSWoPw0UN\nli5lzK3eJUmaOfo+XHT6VNRRK1fWs11Jktqt78NFpy//PWrt2nq2K0lSu/V9uNi7t57teqVOSdJM\n1ffhoi6dvtW7JEmdYrioyfXXe8aIJGlm6vtwUdcVM++4wzNGJEkzU1+Hi8FB2LWrvu1fcIG9F5Kk\nmafvw0Wd9/pYssRwIUmaefo6XEiSpOoZLiRJUqUMFzUbGqq7AkmSqmW4qJmXAZckzTR9HS66YTDl\nVVfVXYEkSdXq63CxeXPdFcCPf9wdIUeSpKr0dbjoBsPDXkxLkjSzGC66gOMuJEkzSV+Hi02b6q6g\ncPPNdVcgSVJ1+jpcdMOYC4DvfhduvBHOOsvxF5Kk3jer7gLq1C0/5KtWwRlnwNKl8Oxnw5w5dVck\nSVLr+rrnYvXquivYZ+nS4q/jLyRJva6vw8X69XVX8GBr19ZdgSRJ09PX4WJ4uO4KHuyii+C007rn\nkI0kSVPV1+Fi9+66K3iw9evhggvgZS8rBnlKktRr+jpc7NlTdwXNDQwUgzzf9jZ7MSRJvaWvw8X2\n7XVXMLGlS+GccwwXkqTe0tfhYmSk7gom56/+ykMkkqTeYbjoAQMDcPXVdVchSdLk9HW46OYxF2Mt\nWeLhEUlSbzBc9IilSw0XkqTe0LfhYmgIdu2qu4qpueWWuiuQJOnA+jZcbNhQdwVT9773eWqqJKn7\n9W246JbbrU/FzTd7aqokqfv1bbjo5Xt4vPrVnpoqSepefRsu7r677gpat2oVnHmmPRiSpO7Ut+Hi\n3nvrrmB6rrwSli+vuwpJkh6sb8NFL465GGvlyrorOLDBQTjrLHtZJKmf9G246OXDIqM+//nO/2hP\nNSwMDsLZZxsuJKmf9G24GBqqu4LpW768uHJnJzWGhbFBw14KSRL0cbjYurXuCqrxiU/U82M+NFQE\nibPP3hcoRoPH/Pnwve91viZJUneYVXcBmp5f/Qp+8zfhj/8YPvIRmDOnM9s9/fTirBWACy6A5z8f\n/uVfitd33w0f/GBxY7hLLulMPZKk7tG3PRczybZtcOml1Z49MjgIixbBz3627/Xb3gbvf3/xejRY\njLrqqsZDNItZuhRe/nJYvLiY8u5377u6qIdPpm/x6I5V27iP28v923517uMZEy4i4skRcU1E3BoR\nP4mIY+quqdNe8QpYsKCaC2wNDBSHXH73d+Hkk+EtbymuDtpsjMf+04sPdOMt7a+8slj+qU+FE04o\nDp8sWVIEDi9pPnV+Mbef+7i93L/tV+c+nkmHRT4HnJ+Zl0TEq4GLgRNqrqmjdu4sQsE55xT3IZk7\nd/LLDg7Cxz4G990H27fD1VcX0/fuhcsuO/Dykx3DsmVL8QC47rqixwWKULRiBZx2WucO7UiS2mNG\nhIuIOAxYALwQIDO/HhGfjognZubt9VbXeZdcUjze/GZ47GPH/8FeswbuuQc+9Sn4/vfhmGOK27p3\n0miwADj3XLj++qI34+Mfh2uuMWhIUq+aEeECOAIYzMy9DdPWAHOBvgsXoz7zmeLvBRfAhRcWQePz\nny/uCHvFFfvPW/d1P66/vvj785/D619f1PiJTxS9Lx//eDFYdPlyOO88OO64Yt7RgHTBBRBRHGpp\nFkbWrCn+jvbmjH3dCXVsU5LqMFPCxVQcXPy5ud4qOmhwEF7ykk5ucRgYaHnpDRvKtQwXYeMP/3Bf\n2wknFGM5Djpo/zEdAF/5CjzykcWhHYDZs+Hgg+HJTy56ZXbtKs5qgX29NHPmwCMeAUccUWz39NPh\npptg3ryiV+cv/gK+853ijJzXvKZ4vm5dEXSe+tTiMvJPe1rR9sMfFutv/PuqVxXbOf98+Pa3i+fn\nnAO/+MW+tiVL9p//sMOKU31Hz7Q55ZR905YsKeZZv36Y004beGAdY+etQuP2mq1zMvM0ztuOOidj\nKnWOGh54+tIIAAAGX0lEQVQeZmCg9c9xu7TyXrrN0BDcdNMwV1010LPv4UCq+Hea7jqq+gw31rFy\n5QO/nQdPtExk5rQ3XLfysMhK4NDR3ouIGAR+f+xhkYh4PfDlzlcpSdKMcXJmNh2RNyN6LjJzKCIG\ngFOAiyPiNcDaJuMtrgJOBlYBOzpXpSRJPe9gYB7Fb2lTM6LnAiAijgYuAn6dol/+jZn5y1qLkiSp\nD82YcCFJkrrDjLmIliRJ6g59FS68iufURcTDI+KKiLglIm6IiKsi4kll22ERcWVErIiI5RFxYsNy\nLbX1s4h4Y0TsjYiXl6/dvxWJiIdFxD+V++SmiPhSOX3sd8LTGpaZSlvff5dExEsiYln5PbE8Iv6s\nnO7nuEUR8cmIuKP8Xji2YXpVn9v2faYzs28ewA+AU8rnrwZ+WndN3f4AHg68uOH1m4Gl5fMvAP+n\nfP4sYC1wUPn6wlba+vUBHAlcUz5ePp196P4dd/+eC3yy4fXh5d+m3wmttvXrA7gXeHr5/EjgfuAR\nfk9Ma58+B/gNius1HdswvfLPbdWf6dp3Xgf/kQ4DNgMPaZg2CDyx7tp66UFxJdTby+dbR7+ky9c/\nAU5qoe260bZ+fAABfA94JrCUfeHC/VvN/v01ikHejxwzvel3Qqttdb/XmvfzEPCc8vmxZRh4qJ/j\nSvbtHaPhoh2f23Z8pmfEqaiT5FU8q/FW4BsRcSgwKzPXN7StAua20Laa4t+hX50B/Dgzb4gIANy/\nlXoSsBH4h4h4AbAdOJviy7TZd8KWFtv6+bvkT4ErImIb8GjgVcAh+Dmu2kS/Za1+biv/TPfVmAtN\nT0ScSfFFfWbdtcwUEfF0ii7I99ddyww2i6Kb/heZeTxFQL68nB51FjZTRMRBwLuBV2bmPOAFwKW4\nj/tWP4WLtcCciGh8z3Mp0pkOICLeDrySYvzFjszcCIxExOENs80DVrfY1q//DidS/PCtjIg7gN8F\nLgBeh/u3KmuAPcBlAJl5I0U385HA45p8J0z0feF3yYMdB8zJzGsAMvN64E6KwyO7/RxXqtXPZkc/\n030TLjJziOKGF6cAxMRX8VSDiDiDosvzhZnZeHP1rwF/U85zPMXAo//XYtuP2vsuulNmnp+Zj8/M\nJ2bmURTHlf8qM8/H/VuJzLyXYrDaiwEi4iiKH6qrafKdMNH3hd8l4xr9cXoqFGceUBzLvwU/x5Vq\n9bPZ8c903QNVOvkAjgb+A7gV+CnlyGYfE+6zxwN7Ke7dMgDcAFxbth1OcQnYFcDPgec2LNdSW78/\ngH9n34BO9291+/Woct8uLz/DryynN/1OaLWtXx/AnzTs35uAPymn+zlufZ+eTxHcdlEMsFxRTq/8\nc1v1Z9ordEqSpEr1zWERSZLUGYYLSZJUKcOFJEmqlOFCkiRVynAhSZIqZbiQJEmVMlxIkqRKGS4k\nSVKlDBeSpi0ivhMRbxpn+o0R8coJlntPRJzT3uokdZrhQlIVLgRObZwQEc8CHgf8ay0VSaqN4UJS\nFb4FPCEifqth2huBLwHHRMSPI+L6iPhFRJw53goi4s8j4oqG1y+NiKUNr98QEdeV6/lhRBxbTv+d\nctpARCyPiNPa8xYlTdasuguQ1PsycyQiLgH+B7AoIh4OLASeTXHr7ZMyc3dEHAz8R0R8PzN/Ot6q\nxnsdEb9Xru/Ecj3PobiF+m8B7wI+mpmXl/PObsNblDQFhgtJVfkC8MOI+F/Aq4H/zMxbI+Iw4LyI\nOI7iDrtPAI6juPPiZL0COBb4SUREOe3RZYhZCvzviDga+PfMvKai9yOpRR4WkVSJzLwZ+BXwcopD\nIv9cNn0AGAKekZnHAT8CDh5nFSPAQQ2vG+cJ4OLMnJ+ZzywfT8jMnZn5SeClwDrgAxHxmUrfmKQp\nM1xIqtIXgDOB44GvltMeA9yZmRkRTwFe2GTZXwHHRsTDI2IW8PqGtm8Bb4iIIwCisKB8fnRmrs7M\nCymCzO9U/q4kTYmHRSRV6XLgXOArmbm9nPY+4JKI+FNgFfCD8RbMzJ9ExHeBnwN3AAOUQSEzr46I\nvweuiIiDgIcB3wGWAW+JiJOAXRS9H2e06b1JmqTIHDt+SpIkqXUeFpEkSZUyXEiSpEoZLiRJUqUM\nF5IkqVKGC0mSVCnDhSRJqpThQpIkVcpwIUmSKmW4kCRJlTJcSJKkShkuJElSpf4/sryeom7yYmkA\nAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1cbc0b991d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# With the series  freqs and vals, a function will plot histograms for bands b5, b3 and b2\n",
    "# # function to plot bands red \n",
    "\n",
    "def plot_hist(vals, freqs, min = 0, max = 6500, title='Histogram', color = 'red'):\n",
    "    \n",
    "    # ind = np.arange(len(vals))  # \n",
    "    #plt.bar(vals,freq, color='#d62728',edgecolor='#d62728')  \n",
    "    #plt.ylabel('Count')\n",
    "    #plt.title('Histogram Band 5')\n",
    "    f, ax1 = plt.subplots(1,1)\n",
    "    #ax1.bar( vals['b5'],freqs['b5'],color='red',edgecolor='red') #d62728\n",
    "    ax1.bar( vals ,freqs,color=color,edgecolor=color) #d62728\n",
    "    ax1.set_title(title, fontsize=9)\n",
    "    ax1.set_xlabel('Values', fontsize=8)\n",
    "    ax1.set_ylabel('Count', fontsize=8)\n",
    "    ax1.tick_params( labelsize =8)\n",
    "    ax1.set_xlim(min,max)\n",
    "    plt.show()\n",
    "    \n",
    "# Call the function    \n",
    "plot_hist( vals['b5'],freqs['b5'], 0, 9114,\"Histogram Red band\")\n",
    "plot_hist(vals['b3'],freqs['b3'], 0,8531, \"Histogram Red band\", color =  \"green\")\n",
    "plot_hist( vals['b2'], freqs['b2'], 0, 10486,\"Histogram Blue band \",  \"blue\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhAAAAF7CAYAAACD27KNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3XucZGV95/HPT0ZlNTpKxIDKcHNx0QRxRmO8xrC6GpOg\nq+LaAiGYGDQm0VHiKtEVojHelokbo8x45WZHXRg1MUgSmRBFRZkGRg0wKJcB6cCMQIO6iDC//eOc\nds60Xd11qqvq1Kn+vF+velXVec6p+p2uvnz7Oc95TmQmkiRJddyn6QIkSVL7GCAkSVJtBghJklSb\nAUKSJNVmgJAkSbUZICRJUm0GCEmSVJsBQpIk1WaAkCRJtRkgpAGLiDdHxNlN1zGqIuJtEXFLRNwR\nEQ9tuh6AiFgZETsjYlWH9o9HxKlDrGddRHxsWO8ndcMAIS1BRGyKiD+dZ/nOiDgMIDP/KjOP7uK1\n3hYR5w6izmGKiGsj4sgu130U8BZgdWY+ODNv6/E9d0bED8sQMh0RkxHxsF5eq8J5/qUFGCCkwRi5\nPz4RsUfTNczjQODOzLxxvsYaNSfwlMx8MPBYYG/gfUusLZa4vTTWDBDSYPzsj0/Zs7Cx8vzd5X/J\nMxFxZUQ8PyJeAJwE/HZE3BkRd5TrroiIv4qI6yPi5rn/WUfEYyPia+Vrfal87U2V9p0R8ZqI+Bbw\nw4h4YESsjYit5X/rV0fEayrr719uc3xEfK9c510RsU9E/FP5Ppsi4uFdfREijouISyPiLWX90xHx\n2rLtBcD5wMryff6ll5orX+8AKHsxzgXWzKnl2RFxcUTcFhHfiojfqbTdLyI+FBE/iIjvAS9h8RC4\nMiI2lp/XZRHxtMrrHV2+xx0RcV1E/MWcWnZGxAnlOrdHxGcj4kGV9mdGxJZy+/8LPAhpxBggpOFI\ngIh4DjABHJ6ZK4FnA1sz83PAO4F/yMwHlf9JQxEqng88leK/dYCzy9daAXwe+ALwi8CbgVfw83/4\nJsr3eTDwY+A64Fnle/wB8N6IeMqcbZ4FPA74VeB1wKeAPwUeBvy0rKtbjwN+CDwCeFn5fgeW+/yb\nwEx5+OLZS6wZgIj4JeAo4KrKssOATwNvzMyHAicAZ0bEfy5XeQvwZIreiycAL+pivyaADwMrgQ8B\nn4+I2c9tB/Dfy3qPBF4ZERNztj+K4uu8CtgPWFvW+hDgc8D/AR4CfAI4pot6pKEyQEhL966IuLVy\nu43O/73+FLg/8CsRsSIzb8zM7y7w2scAb8/M72fmj4HXA8+JiH2ApwB7Ae/MzHsy8xsUf+jnendm\n3pyZP83Cxsy8CSAzL6ToBXjWnG3enpl3ZeaVwOXAVzLzysz8KbARWN3F12XW9sz868y8t3y/64DD\nF9mml5q/XPbcTAN7An9SaftD4OPltmTmV4F/AF5atr8c+MvyPe8ATulivy7IzH/MzJ2ZuR64Gfjt\n8vXPn/1cM3ML8Hfz1PvuzPxB+X7nsKvH5HeA72fmR8rX/gfggi7qkYbKACEt3Zsyc6/K7aF0OH6e\nmf8KvA14O7A9Ij4TEfsv8NqPAq6vbD8N/KRcvi8wnZk7K+tvm+c1bqg+KbvXN5fd9bdR9AJUBxwm\ncEvl+Y8p/jhWn//CAjXPdfOc5z9i8S75ujUDPL38j//XKL4+j6i0HQC8ak7IO5Lia0i5bvVrdz2L\nj4G4fp7njyzrfW5EXBQR2yPidooej7n1Vr8u1a/Jvh1eWxopBghpyDLztMx8CkXX9d0UXdUAO+dZ\n/UaKP34AlD0P9yuX3wTsExHVn+P5Tjv82etGxH4UXeInAg8rw855jN6AwV5qnh0D8Q3gfwMbKm03\nAO+vhrzysMkfl+03AdUgtz+Lj4GYG/xWATdGxH0pehQ+BOybmQ8B1s9Tbydza5l9bWmkGCCkIYqI\nJ0bEU8o/Mj+h+M/znrL5ZmD/2P3Mg7OAkyLiURHxC8CpwD9n5n8AXwduB95cDrZ8Eru65DuZ7TnY\nXtbzfOC/zS2zx90blG5qnus04FERMTuWYT1wfEQ8KyLuExH3j4hfi4jHlO2TwJsiYt9yDMJbu6jr\niIj4zYjYIyJeCewD/CPFIar7A7dm5j0R8WSKQyTd+gLwyIj4/fK1fws4osb20lAYIKSl6fRfaqfl\nDwY+SDHI7iaK7urXlm2fAe6gOLRxa7nsryiO938NuAbYAzgWIDPvAV5Accz8VuBdwJkUwWTeOjLz\nCuAvgU1lDUdRDNhbqPa6p6Qutn6t9l5qzsy7gHXAyeXzyygGPb6DIojcAPwFxR96yuWXAN8GpijG\neSxW4yeBV1KEuD8GjszMmcz8IfAa4MPl4Ys3U4yB6FjvnNpvo/hcXwfcRjEw9qxF6pGGLjJH7nR1\nST2KiNMofq5PaLoWSePNHgipxSLi6eXhjYiI/0rRVf7ppuuSNP5WDPsNI2Iv4Evs6sJ7IMX57Q8H\n7gucARwM3AW8JjO/XG63d6c2aRk7iKJ7/CEUAyv/Z2Z+qdmSJC0HjR/CiIg3AM/MzBdEcbGY6zLz\nLyLiiRTHIQ/IzHsj4qPA9fO1NVi+JEnL0igcwvh94CPl46MoRk+TmZdQDDL79bLtpXPavl9pkyRJ\nQ9RogIiIp1J0vX6hPLSxIjOrE9hcB6zq0HY9nhstSVIjhj4GYo5XAGdk5s6I/p16HhG/CDyXIoDc\n1bcXliRp/O1JMYHd+Zn5g04rNRYgIuKBFIclngiQmbdGxD0R8fBKT8MBFOMeOrXNN20vFOHh7MFV\nL0nS2DuaYr6TeTXZA/Ey4LLM3FpZ9hng1cAp5ax6jwD+bYG2Czu89nUAZ511FoceeugASm/G2rVr\nWbduXdNl9JX71A7uUzu4T+0w6vt0xRVXcMwxx0D5t7STJgPE8ew+Vz3AmygusbuVYja9oytnWSzU\nNtddAIceeiirV9e5aOBoW7ly5VjtD7hPbeE+tYP71A4t2qcFhwA0FiAy8+nzLLuF4vDDfOt3bJMk\nScM1CqdxSpKkljFASJKk2gwQLTIxMdF0CX3nPrWD+9QO7lM7jMs+NT6V9SBExGpg8+bNm9syUEWS\npJEwNTXFmjVrANZk5lSn9eyBkCRJtRkgJElSbQYISZJUmwFCkiTVZoCQJEm1GSAkSVJtBghJklSb\nAUJqyDRwcnkvSW1jgJAaMg2cggFCUjsZICRJUm0GCKkh2+fcS1KbGCCkhuyYcy9JbWKAkCRJtRkg\npAZMA+c2XYQkLYEBQmqAAUJS2xkgJElSbQYISZJUmwFCkiTVZoCQJEm1GSAkSVJtBghJklSbAUKS\nJNVmgJAkSbUZICRJUm0GCKlh51LMTClJbWKAkBpmgJDURgYIqQHbmy5AkpbIACE1YEfTBUjSEhkg\nJElSbQYISZJUmwFCkiTV1kiAiIj7RcTfRMTWiLg8Is4olz86Ii6KiKsi4uKIeGxlm7lthzZRuyRJ\naq4H4t3Azsw8JDMfD5xYLl8PnJaZjwHeA3yiss3cttOHWK8kSaoYeoCIiAcArwD+fHZZZt4SEXsD\na4Czy2XnAPtFxEELtQ27fkmS1EwPxMHArcCfR8Q3I+LCiDgC2A+YzsydlXW3AasWaZNaZRo4p+ki\nJGmJmggQK4D9gW9n5pOA1wKfKpdHA/VIQzUNbGy6CElaohUNvOc24F7gkwCZeVlEXEsRKvaJiPtU\nehpWlevfCezboa2jtWvXsnLlyt2WTUxMMDEx0bedkSSprSYnJ5mcnNxt2czMTFfbDj1AZOYPIuJL\nwPOA8yLiQOAA4CvAFHAscHpEvAS4ITOvAYiIjm2drFu3jtWrVw9uZyRJarH5/qmemppizZo1i27b\nRA8EwKuBj0bEuyl6I/4wM6cj4lXAJyLiJGAGOL6yzUJtkiRpiBoJEJl5LXDEPMu3Ak/tsE3HNkmS\nNFzORClJkmozQEiSpNoMEFKLTQMnl/eSNEwGCKnFpoFTMEBIGj4DhDQG1mOIkDRcBghpDGzAACFp\nuAwQkiSpNgOENAI8BCGpbQwQ0gjwEISktjFASC22vekCJC1bBgippaYpei4kqQkGCKmlpoFzmy5C\n0rJlgJAkSbUZICRJUm0GCGnI+nW2hQMoJTXJACEN0TbgRX16rR19eh1J6oUBQhqiHcDdA3pteyQk\nDZMBQhoTVzZdgKRlxQAhjYk7my5A0rJigJDGhAFC0jAZIKQRNg2cTHdnblzQ5XqS1A8GCGmIFhvo\nOA28obxNl7dT6C4YXNLlepLUDyuaLkBaThY79XIaOLV8fPSAa5GkpbAHQmohr4MhqWkGCGlErGfh\nQxzV8RAGCElNM0BII2IDCx/iqDMeQpIGzTEQ0ojazq5ehupjSRoF9kBII2oHRa/E3McLWY89FJKG\nwwAhjZENGCAkDYeHMKSW8XCGpFFgD4TUMlfT3eEMSRokA4Q0om7rsPz2oVYhSfMzQEgjyqAgaZQZ\nIKQh6tSr0G17Nxa73oYk9YMBQhqixXoV+tHrsNj1NiSpHwwQkiSptkYCRERcFxFXRMSlETEVEUeV\nyx8dERdFxFURcXFEPLayzdy2Q5uoXWrCOZXHdzZWhSTt0lQPxE7gpZn5hMxcnZmfKZevB07LzMcA\n7wE+UdlmbtvpwyxYGrZqUNjYYbkkNaWpABHlbdeCiL2BNcDZAJl5DrBfRBy0UNtQq5aGyKAgaZQ1\nOQbijIi4PCI+HBEPA/YDpjNzZ2WdbcCqRdqksXFBF+tcNPAqJGlxTU1l/YzMvDEi9gD+kuJwxFuZ\n0yuxVGvXrmXlypW7LZuYmGBiYqKfbyP1zaYu1tmySPu5wBHAvksvR9KYm5ycZHJycrdlMzMzXW0b\nmTmImroWEfsAVwEHA98DHjrb0xAR08DTKHpzrwb2mtuWmdfM85qrgc2bN29m9erVw9kRqQvvoEjK\n3Xg18KEe32cz4He+pF5MTU2xZs0agDWZOdVpvaEfwoiIB0REtVvg5cBUZu6g+L13bLneS4AbMvOa\nzNwOTM3XNtzqJUkSNHMI45eAcyLiPhSHLK4BjivbXgV8IiJOAmaA4yvbLdQmSZKGaOgBIjOvpUPv\namZuBZ5at02SJA2XM1FKkqTaDBCSJKk2A4QkSarNACFJkmozQEhDMk13M032w/YhvY+k5csAIQ3J\nNN3NNNkPO4b0PpKWLwOENKK+2XQBkrQAA4Q0oi5pugBJWoABQhqCaWB900VIUh8ZIKQhmAY2NF2E\nJPWRAUKSJNVmgJDGkGdhSBo0A4Q0hk4EtjVdhKSxZoCQxtA92AshabAMEJIkqTYDhCRJqs0AIUmS\najNASJKk2gwQkiSpNgOEJEmqzQAhjantTRcgaawZIKQx5TwQkgbJACFJkmozQEiSpNoMEJIkqTYD\nhCRJqs0AIUmSajNASEPgKZWSxo0BQhoCT6mUNG4MEJIkqTYDhCRJqs0AIUmSajNASJKk2gwQkiSp\nNgOEJEmqrdEAERHHR8TOiDiyfL53RJwXEVsjYktEPKOybsc2ST/vHGC66SIkja3GAkRE7A/8AfC1\nyuJ3AV/LzEOAVwCfjIg9umiTNMdGDBCSBqeRABERAXwE+GPg7krTS4HTADLzEuAm4Nc7tH2/0iZJ\nkoaoqR6I1wNfzsxLZxdExF7Aisy8pbLedcCqDm3XA6uGUay0FNPAuU0XIUl9tmLYbxgRjwNeDAx8\nDMPatWtZuXLlbssmJiaYmJgY9FtLP2OAkDSqJicnmZyc3G3ZzMxMV9sOPUBQBIf9gavLQxn7ABuA\nk4F7IuLhlZ6GA4DrM/PWiJivbdtCb7Ru3TpWr149gF2QJKn95vunempqijVr1iy67dAPYWTmaZn5\nyMw8KDMPBL4OvDIzTwM+A7waICKeBDwC+Ldy0/naLhx2/ZIkqZkeiLkSiPLxm4AzI2Ir8BPg6My8\nt4s2aWR5KW9J46jxAJGZR1Qe3wI8t8N6HdukUealvCWNI2eilCRJtRkgJElSbQYISZJUmwFCkiTV\nVjtARMRnu1kmSZLGVy89EPNNH33QUguRJEnt0fVpnBFxAvAq4JCImKo0rQS+0+/CJEnS6KozD8QX\ngauADwFrK8vvALb0syhJkjTaug4QmXk9xRUwDx1cOZL6aT3FRWb2bbgOSeOn9kyUEXEA8D+Bg6vb\nV2eUlDQaNgAnYICQ1H+9TGX9aeBLwAcAr0UhSdIy1EuA2DMz39z3SiRJUmv0chrntyNivlM5JUnS\nMtFLD8TewOUR8TXgrtmFmfmivlUlqe+mKQZVOiZCUj/0EiDOKm+SWmQaOAU4EgOEpKWrHSAy8/RB\nFCJpMGZP5ZSkfurlNM6Pzbc8M1+x9HIk9dvsqZyS1E+9HMLYXHm8J/BiYKrDupIkaQz1cgjjb6vP\nI+JDwGf6VpEkSRp5vZzGOdddwIF9eB1JA7K96QIkjZ1exkCcWnm6B/BE4N/7VpE0RqaBs5suArga\nOLfpIiSNlV56IGYqt1uA/wNM9LMoadRNU5zZMN3FOucNoZ7F3E4xmFKS+qWXMRCnDKIQqU26mVNh\nGv9oSxpftXsgIuJBEfG3EbG1vH0gIh40iOIkSdJo6uUQxgcpei5eChxVvsYH+1mUJEkabb3MA3FY\nZj6+8vyPIuLyfhUkqf8uaLoASWOnlx6IPaqHLMrHe/SvJEn9tqnpAiSNnV56IE4Hvh4RnyqfvxT4\neP9KkiRJo67rHoiIeHBEHJCZ7wVOBB5c3j4IfHhA9Unqs/UsfPqpJHWjziGM9wBrADLzvMw8MTNP\nBO4E3j2I4iT13wYWn8NCkhZTJ0A8MzPPmbswM88Entm/kiQN2gYMEJKWpk6AuHeBtp1LLUSSJLVH\nnQCxIiIePHdhRKwE7tu/kiRJ0qirEyD+DjgzIh46u6B8/PGyTZIkLRN1AsQ7KK7Jc0NEXBoRlwI3\nUAyifPsgipPaopuLa0nSOOk6QGTmvZl5HPB44J3l7fGZeVxmLjQ+4udExPkRcVkZRC6MiMPL5Y+O\niIsi4qqIuDgiHlvZZm7boXXeUxqk2YtrGSAkLRe9XI3ze8D3lvi+R2XmHQAR8ULgE8DhFKeon5aZ\nZ0bEi8vlv1puM7ft9EqbJEkaol6msl6y2fBQegiwMyL2Bp4InF2ucw6wX0QcVLatma9tuJVLu9tO\ncehie8N19KKNNUsaHb1MZd0XEXE68BtAAs8H9gNuyszqKaHbgFXAHcB0h7ZrhlOx9PN2UBy6OKvp\nQnqwo+kCJLVaIz0QAOXYiVXAWyhmuQSIpuqRJEnda6wHYlY5pmE9xRkd+0bEfSo9DasoehruXKCt\no7Vr17Jy5crdlk1MTDAxMdHXfZAkqY0mJyeZnJzcbdnMzExX2w49QJQTTz0gM6fL5y8EdmTm9ojY\nDBwLnB4RLwFuyMxryvWmOrV1sm7dOlavXj3I3dGYmaYYrXsCsG/DtQzLctxnSYX5/qmemppizZo1\ni27bxCGMlcBnI+LyiLgM+CPgt8u2VwEnRMRVwBuB4yvbLdQm9UW3p2PODkA8e87y2StdTgPv629p\nA+MpqJJ6MfQeiMzcBjy5Q9tW4Kl126R+W09xdkWn/8hnByCeN2f5Bor/5KeBSSRpfDU2iFIaZUu5\nWuV64No+1jJIs4cvJKkuA4RU0zRw7gLtGyjmfB915wBbKOqVpLoMEFJNiwWIttiIc0FI6p0BQupg\nKTM1XtC3KgbrnKYLkNRaBgipg6X8d76pb1UM1samC5DUWgYIqQ9ua7oASRoyA4TUpWk6XzjrhuGW\nIkmNM0BIHcztVZidcGm+Qxt3Dr4cSRopBgiptI3d536ocyrmN/tciySNusYvpiWNgm3AY4Dq9eI/\nSjHH+uFdbH/JIIqSpBFmD4REcVjiLuDuyrLrgO80Uo0kjT4DhFTTOJ9xMTtQ1AtrSVqMAUKiu0mj\nZtcZ54tkeWVOSd0yQEh0N2nUleX9VwdZSIPWs7TZNyUtLwYIqQvbgD9ruogB24DXxpDUPQOE1IUd\nwE+bLkKSRogBQurCcunaP7vpAiS1hgFC6sJy6do/r+kCJLWGAUJaxGXAO5ouQpJGjDNRSgs4l2Le\nhysXW1GSlhl7IKQFnFvelpv1OBeEpIUZIKRFbGq6gAZswAAhaWEGCEmSVJsBQpIk1WaAkDSv5TL3\nhaTeGCAkzWu5zH0hqTcGCEmSVJsBQpIk1WaAkDSvc/FUTkmdGSAkzetcYEvTRUgaWQYIiWK6av08\nB1JK6sQAIQG3N12AJLWMAUKSJNVmgJDUkYd2JHVigJDUkYd2JHUy9AAREfePiI0RcWVEXBoR50fE\nwWXb3hFxXkRsjYgtEfGMynYd2yRJ0nA11QOxPjP/S2Y+Afg88JFy+buBr2XmIcArgE9GxB5l27sW\naJM0AHc2XYCkkTX0AJGZP8nML1YWfR3Yv3x8FHBaud4lwE3Ar5dtL53T9v1Km6QB+CLwBpxQStLP\nW9F0AcBrgc9GxF7Aisy8pdJ2HbCqQ9v1wKrhlSktP1vK2w+Bk4F9G61G0ihpdBBlRJwEHAyc1GQd\nkha2AXshJO2usR6IiDgReCHwXzPzLuCuiLgnIh5e6Wk4ALg+M2/t0LZtofdYu3YtK1eu3G3ZxMQE\nExMT/dwVSZJaaXJyksnJyd2WzczMdLVtZOYgalr4TSNeD7ycIjzMVJZ/jCIwnBIRT6KYjv+AzLx3\nobZ5Xn81sHnz5s2sXr16KPukdnsH8NamixhxmwF/mqTxNzU1xZo1awDWZOZUp/WG3gMREY8E3gd8\nD9gUEQHclZlPAd4EnBkRW4GfAEdXAsJCbZIkaYiGHiAy8/t0GHtRHp54bt02qVfTwHrg/zVdiCS1\njDNRalmbBk7B+Q4kqS4DhAT8qOkCJKllDBAS8OOmC2iJaYr5IDylU5IBQgK+03QBLTF7yMcAIckA\nIQFXNF1AC6wDtjddhKSRYYCQ1JWzgB1NFyFpZBggJHXttqYLkDQyDBCSunZ70wVIGhkGCEmSVJsB\nQsuagwIlqTcGCC1rDgqUpN4YILSsOSiwntkpv+25kWSA0LI1DUw2XUTLXFPeX91oFZJGgQFCy9Ls\nlMxfbbiOtvmX8t6zMSQZILTszIaHDQ3X0UYGB0mzDBBadqYxPEjSUhkgJNVmT4QkA4Sk2v4a2NZ0\nEZIaZYCQVNu9OIeGtNwZICRJUm0GCEmSVJsBQpIk1WaAkNQTp7OWljcDhKSevJtiTg1Jy5MBQlJP\nNgFb2DWzp2FCWl4MEFpWpoH1TRcxRq6k+JqeggFCWm4MEFpWtuA01v30ORwLIS1XBggtK05+1F+b\n2HVp71OxF0JaTgwQkpbkhvL+bIoeHknLgwFC0pLcWXlsD4+0fBggtGxMA+c0XcQYurnpAiQ1wgCh\nZWMa2Nh0EWPolqYLkNQIA4SkJflh0wVIaoQBQtKS3Nh0AZIaYYCQtCQOnJSWJwOElg0nPBq8s4HL\ncGpraTkYeoCIiPdHxLURsTMiDqssf3REXBQRV0XExRHx2AXaDh123Wo//1MevPOA1+PU1tJy0EQP\nxGeApwHXzVm+HjgtMx8DvAf4xAJtpw++TI2b25ouYJnY1HQBkoZi6AEiM7+SmTcBMbssIvYG1lD0\ngJKZ5wD7RcRBC7UNu3a12+1NFyBJY2RUxkDsB0xn5s7Ksm3AqkXapEXNXm76zkXWU39d2XQBkgZq\nRdMFDNLatWtZuXLlbssmJiaYmJhoqCI1YQvFMfnDFltRffW3wG8A+zZdiKSOJicnmZyc3G3ZzMxM\nV9tGZg6ipsXfOOJa4AWZuaU8THE1sNdsT0NETFOMlbizU1tmXtPhtVcDmzdv3szq1auHsDcaZWcD\nxzRdxDL1IuADGCKkNpmammLNmjUAazJzqtN6I3EIIzO3A1PAsQAR8RLghsy8ZqG2purVaJs9ZOFZ\nAM07FwdVSuOqidM4T4uIG4BHAudHxNay6VXACRFxFfBG4PjKZgu1SbuZZvfTCD37olmXN12ApIEY\n+hiIzHxVh+VbgafWbZM6WU/RE+HZF836j6YLkDQQYz2IUsvbhvLeb/Jm/aDpAiQNxEiMgZAGZQPw\nT00Xscx9FXgKxRTXjk+RxocBQmPvu00XsMzdBnwd+Dw/Pz5FUnsZICQNxV+wKzisxxAhtZ0BQtJQ\n3Au8rny8AQOE1HYGCI2N2ePrXrZ7dFUPJ/k5Se1mgNDYmD2+7mW728HPSWo3A4SkRlxX3ntmhtRO\nBgiNnXObLkBd+TRFaJjtOfoz4A0YJKS2cI4djR0DRDtsoQgNUT4/u7w/Gi++JbWBAUJSY85efBVJ\nI8pDGBoL0xRzC0iShsMAobEwza5rX6jdLgTW4lgIadR5CEPSSHkDkMBq4NiGa5HUmT0QkkZKlvcf\np7gA18nYGyGNIgOEWm12DoGrGq5D/bcJeDtefEsaVQYItdrsHAIfbroQDcTsKbnbsSdCGjUGCLXa\n7PUUNjVahQZtB/ZESKPGAKFW83oKy8N15b2XAZdGhwFCrXZb0wVoKN5b3m/AQxnSqDBAqNVub7oA\nDcVM5fEG4AiKMzQkNccAodaaBj7XdBFqxJXAd5ouQlrmDBBqrU3AJU0XocZ8FK/eKTXJmSg18mav\nc3FC+fx9wI+AbzVWkUbBpvJ2M8UYieoVPKvfM17ZUxoMA4RG3uxcD98tn3sFR1WdTXH9jF8FDgBO\nZNf3zJEYIKRBMUCoNQwO6uTG8gZwBfCLDdYiLRcGCI2MuYcq1gNPA97UWEVqo/Mqj7d3XEvSUhkg\nNBKmgT+j6GX4FvCfyse/C0w1WJfa7WLguXOWOT5C6g/PwtBAzF7kanqeZdUrLM4u+1d2HaI4t/L4\njAHXqfH2NuBF7H62xuz4iOrzk/FsDqkueyA0EPMNYptd9rDy/nvAfysfP6qBGrU8bCzvL6WYkOqX\nyudvoTiD4wCK0OqAS6keA4SGbvYKi2cBnywf39hhXalf5l5wbXasxOwhMsdLSPUYILRks8eUXwh8\nll2DIAF+D3g88AB2Xfiq+ot85xDqk7rxduCw8rFjJKTFGSDUk20U3b9fAB7D7ocmPsWuaxd8Cyd8\nUjtcRBF2H0cxJucQiutu/DVweHNlSSPLAKGfmTs6/TLgdcCfA18s13ke8FaK48kAdwMPKh9/qLy/\nchjFSgP6ng88AAAKa0lEQVSwnSI8ALyWotfsd4DfohiQ6WyX0i6ehdGwaYoR4ovN6b8N+JvJyb69\n7z9TDB77ALCWXWdEnFJ5fBLFDH+fBU4tb8dRnBp3d3kDuLO8//deCunjPo0M96kdFtmn2UNuN1IE\nhWOA1cDpFD+Pc8/mGAWTY/g5uU+jq1UBIiIeHREXRcRVEXFxRBzadE1LNc2uP84nU/zXfwLFL6rZ\nQV7bKA4TvG5ykm2LvNbJ7AoAs8Hkm8Dfs3tQOBm4HvgTii7aV5TrAhwLPLLy/qfPeY++GpMfpN24\nT+1Qc58uoOh5+z2K8P3ycvlxdL60+LbyNizj8oepyn0aXW07hLEeOC0zz4yIF1P8bfvVhmvqWvWQ\nwD9SnML42Ur7BuCh5T3A84GDgQDuKpcdQTEg8Q8pzmJ4O8Xhhf8Avg1sAf4GuJ1dAxRPrbzHX3eo\n7YuVx3N7En606J5Jy0sCV5WPvw08gWL8xE7gP1OE81XA58t1D6W44NdvDr1SaXBaEyAiYm9gDfAc\ngMw8JyI+EBEHZeY1zVa3uG8Cb6Q4JHAJxR/l+f6Ynzrn+fc6PP+T8v4odg1YnHVr72VK6tHl5f3s\noOHNlbbvUPxD8Bjg14E/KJd/gfnHUMz2WqzqcrnUT92e0tyaAAHsB0xnZvXMv20UP0sjESC+STGv\nwXEUXSMPpRhj8DCKC/zMWug/+p/WfM+54UHS6LqqvG2oLHsvcF+Kwcg3U/Ri3Fu2PQT4DeCVwCco\neizvU77GfCHiduD3gVdRTM42O8hz7u+VVeweRtoWTNpWb9Pqfr12LL4K0K4AUceeAFdcccVi681r\nO8VkRy8C9qb4Yf0rih/I/wT8P+Ba4MDy+c3A19h1yGBuz0LfJqiZmYGpMbsyhPvUDu7TwPy4vJ/v\nn4HbKWbS3Dhn+f6dXmxmhi9PTfGxyqJT5lntgew6LLo3xSFQgD8Cvg68BJhk9995VwP3UPxy/ROK\nHpdnAf+3bJsoH59I0dMy19zfq1D8bn0fRei5vNJWXXdmZobzp6Z+9vwe4GXl9h8EvkJx+OijC7x3\nN+arb1BmZmaYGtL33jS7vl5/R3dnDF2762/nngutF5nZe2VDVB7CuBrYa7YXIiKmgafNPYQRES/H\nqz9LkrQUR2fmJzs1tqYHIjO3R8QUxUkCp0fES4AbOox/OB84GriOXUFbkiQtbk+Kk43OX2il1vRA\nAETEIRSHAn+Rosfv+Mz8TqNFSZK0DLUqQEiSpNHQqomkJEnSaBi7ANGW2Soj4v0RcW1E7IyIwyrL\n59b/2B7bhr7fEXH/iNgYEVdGxKURcX5EHFy27R0R50XE1ojYEhHPqGzXU9sQ9+v8iLis3KcLI+Lw\ncnlrP6uyjuPL778jy+et/YzKOq6LiCvKz2kqIo4ql7f2c4qI+0XE35Rf28sj4ow271NE7FX5fKbK\nOu6OiIe0+fsvIp4fEZvLfdsSEb+7lLpHYZ+6kpljdQO+BBxbPn4x8I2ma+pQ59OBR1DMYXFYN/X3\n2jbEfbo/8LzK89cAm8rHHwP+V/n4icANwB7l84/20jbE/Xpw5fELgcvG4LPan+IClBcBRy7lcxiF\nz6h872uAX5lneZs/p3XA+yvPH972fZqzf28APlc+bvPviB8Ajysf709x5usD27xPXe130wX0+UPc\nm+LU6ftUlk0DBzVd2wI1X0sZIBaqv9e2hvdtDXBN+fjO2V9+5fOLgSN6aPv6bFtD+/R7wFT5NZ9p\n42dFMTv6P1PMwLyJXQGi1Z9R9Wepsqy1P1MUs9bPAL8wLvs0zz7+O/A7bf/+o5hG4unl48Mo/uDf\nt8371M2tNadxdmnkZ6tcxEL139FjW5P7/VrgsxGxF7AiM2+ptF0HrOqhbfYyA0MVEadTTAqYFLMS\n7wfc1NLP6vXAlzPz0ogAiq5lWv4Zlc6IYqe+AbyZdv9MHUwxM/2fR8SzKeacOoUiBLR1n34mIp5K\nMdnmF8bg++9lwMaI+BHFPr2IYnLRNu/TosZuDIRGQ0ScRPEL8KSma+mHzDwuM1cBbwHeUy6OBkvq\nSUQ8jqLr+i+brmUAnpGZh1NczPYH7LqQbOs+p9IKiu7wb2fmkygC+afK5W3dp6pXAGfMCTStExF7\nUPxeeGFmHgA8m+Jah+PyOXU0bgHiBmDfiKjuV3Wq91G3UP29tg1dRJxIMVbgeZl5V2beCtwTEQ+v\nrHYAcH2PbY19npl5JkVPRFs/q2dQ/FG6OiKuBX6N4tIML6Xln1Fm3lje30sxo/zTyzra+DlRvte9\nFJfYITMvozhMsz+wT0v3CYCIeCDF99zHAFr+O+JwYN/MvAggMy8BbqQ4lPHTlu5TV8YqQGTmdorj\n08cCxMKzVY6chervtW3Y+xARr6fozntOZt5ZafoM8OpynSdRDCD9tx7bLhzsXuwSESsjYt/K8xcC\nO8qv+WZa9lll5mmZ+cjMPCgzD6Q4tvrKzDyNln5G5fs+ICJWVha9HJjKzB208HMCyMwfUAx6fF5Z\nw4EUf0i+0kvdo7BPFS+jGIy8tbKsrd9/s8Hsv5Q1PJpizMmVPdQ9KvvUnaYHYfT7BhwCfJXiOi3f\noBwZO2o34DSKb7y7KQYybV2s/l7bhrhPj6S4ptjVFL+oLgW+VrY9nGJa1K0UVzx+ZmW7ntqGtE+r\nKAY3XQ5cBvwTuwa9tvazqtRyAbsGUbbyMyprOLD8nrus/Kw2Aqva/jmV+3UBsKX8eXph2/eprOMr\nwO/OWdbm77//UfmMLgf+R9v3qZubM1FKkqTaxuoQhiRJGg4DhCRJqs0AIUmSajNASJKk2gwQkiSp\nNgOEJEmqzQAhSZJqM0BIkqTaDBCSuhYRX4iIP5pn+WXlNN+dtntbRJw62OokDZMBQlIdH6W4iuLP\nRMQTgX2Av2+kIkmNMEBIquPzwKMi4pcry44HzgAOjYgvR8QlEfHt8pLuPycijouIjZXnvxURmyrP\nj4mIr5ev868RcVi5/MnlsqmI2BIRJwxmFyV1Y0XTBUhqj8y8JyLOBH4fWBsR9wcmgKdQXML4iMz8\naUTsCXw1Iv4lM78x30vN9zwinlq+3jPK13k6xeWsfxl4M/DezPxUue5KJDXGACGpro8B/xoRfwa8\nGPj3zLwqIvYGPhQRh1NclfVRwOEUV33s1guAw4CLIyLKZQ8pg8om4K0RcQhwQWZe1Kf9kdQDD2FI\nqiUzrwC+CxxJcfjiI2XTO4HtwOMz83DgQmDPeV7iHmCPyvPqOgGcnpmrM/MJ5e1RmfmTzHw/8FvA\nTcA7I+Jv+7pjkmoxQEjqxceAk4AnAZ8ulz0UuDEzMyIeAzynw7bfBQ6LiPtHxArg5ZW2zwPHRMR+\nAFFYUz4+JDOvz8yPUoSVJ/d9ryR1zUMYknrxKWAd8HeZ+eNy2TuAMyPiZcB1wJfm2zAzL46IfwS+\nBVwLTFGGgcz8SkS8EdgYEXsA9wO+AGwG/jgijgDupujFeP2A9k1SFyJz7lgmSZKkhXkIQ5Ik1WaA\nkCRJtRkgJElSbQYISZJUmwFCkiTVZoCQJEm1GSAkSVJtBghJklSbAUKSJNVmgJAkSbUZICRJUm3/\nH26dl00jl7MiAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1cbc41fefd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# a plot of infrared band\n",
    "plot_hist( vals['b7'], freqs['b7'], 0, 8766, \"Histogram InfraRed band\",\"cyan\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "It is common to apply a linear stretch, deleting  certain percent of pixels considered noise or outlier. The following example deletes 5% of the values by applying *percentile* reducers to get the percentile 2,5 and 97,5% respectively. These values are  set as *minimum* and *maximum* to scale image values to 0 and 255 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "done!\n"
     ]
    }
   ],
   "source": [
    "# first define a function to transform the image values \n",
    "#=================================\n",
    "# function to transform the image to 255 bits (uint8):\n",
    "def transformTouint8Reducers(image, perc):\n",
    "    NameBands = image.bandNames()\n",
    "    minValues = image.reduceRegion(\n",
    "          reducer=ee.Reducer.percentile([perc]),\n",
    "          bestEffort=True,\n",
    "          maxPixels=100000000,\n",
    "          ).toImage(NameBands)\n",
    "\n",
    "    maxValues = image.reduceRegion(\n",
    "          reducer=ee.Reducer.percentile([100-perc]),\n",
    "          bestEffort=True,\n",
    "          maxPixels=100000000,\n",
    "          ).toImage(NameBands)\n",
    "\n",
    "    return image.subtract(minValues).divide(maxValues.subtract(minValues)).multiply(255).byte()    # this scale values to 0-255\n",
    "\n",
    "# Compute the image with pixel values between 0-255\n",
    "img256 = transformTouint8Reducers(image,2.5)\n",
    "print(\"done!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Compute new histograms\n",
    "geom = img256.geometry()\n",
    "scale = img256.projection().nominalScale()\n",
    "# compute the histogram\n",
    "bfreq = img256.reduceRegion(reducer = ee.Reducer.frequencyHistogram(),geometry = geom, scale = scale,\n",
    "                             maxPixels = 10000000 ).getInfo()\n",
    "\n",
    "nfreqs= dict()\n",
    "nvals = dict()\n",
    "\n",
    "for l in bfreq.keys():\n",
    "    d = bfreq[l]\n",
    "    freq = []\n",
    "    val = []\n",
    "    for k,v in d.items():\n",
    "        if not(k == 'null'):\n",
    "            val.append(int(k))\n",
    "            freq.append(int(v))\n",
    "    nfreqs[l] = freq \n",
    "    nvals[l] = val\n",
    "    #print('Minimum and maximum value for : %s is %.2f %.2f' % (l, np.min(val), np.max(val) ))\n",
    "  \n",
    "# Call the function    \n",
    "plot_hist(nvals['b5'], nfreqs['b5'], 0, 255,\"New histogram Red band\")\n",
    "plot_hist(nvals['b3'], nfreqs['b3'], 0, 255, \"New Histogram Red band\", color =  \"green\")\n",
    "plot_hist(nvals['b2'], nfreqs['b2'], 0, 255,\"New Histogram Blue band \",  \"blue\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Show the image before and after the histogram stretch\n",
    "from IPython.display import display\n",
    "# image before\n",
    "img1  = Image(url = image.getThumbUrl({'min': 0, 'max': 2048,'bands': 'b5,b3,b2'}))\n",
    "# image after\n",
    "img2  = Image(url = img256.getThumbUrl({'min': 0, 'max': 256,'bands': 'b5,b3,b2'}))\n",
    "\n",
    "display(img1, img2)\n",
    "# you may notice that changes are saved into the image img256"
   ]
  },
 
  {
   "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}