earth-observation-for-journalism/true-color-image/02b Timeseries.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "from multiprocessing import Pool\n",
    "import geopandas as gpd\n",
    "from sentinel_helpers import scihub_band_paths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "base_path = Path('input/tempelhofer_feld/')\n",
    "input_files = list(base_path.glob('*.zip'))\n",
    "area_of_interest = gpd.read_file(base_path / 'tempelhofer_feld.geojson')\n",
    "ndvi_path = Path('output') / 'ndvi'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio as r\n",
    "import rasterio.mask\n",
    "import rasterio.plot as rplt\n",
    "import numpy as np\n",
    "from sentinel_helpers import scihub_normalize_range, scihub_cloud_mask\n",
    "from zipfile import BadZipFile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# uncomment this cell to remove previous calculations\n",
    "# ! rm -rf {ndvi_path / '*.tif'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_ndvi(product_path):\n",
    "    global area_of_interest\n",
    "    \n",
    "    try:\n",
    "        b04_path, b08_path = scihub_band_paths(product_path, ['B04', 'B08'], '10m')\n",
    "        cloud_mask = scihub_cloud_mask(product_path)\n",
    "    except BadZipFile:\n",
    "        # This is due to a problem of Scihub\n",
    "        print(f'{product_path}: Problem reading product, skipping it')\n",
    "        return\n",
    "\n",
    "    with r.open(b04_path, 'r') as b04, r.open(b08_path, 'r') as b08:\n",
    "        # we want to only write the bare minimum data necessary to disk\n",
    "        out_meta = b04.meta.copy()\n",
    "\n",
    "        # we reproject the geojson file we fetched above and convert it so that rasterio\n",
    "        # can use it as a mask; we subtract the cloud mask to avoid irregularities\n",
    "        reprojected = area_of_interest.to_crs(out_meta['crs']).iloc[0].geometry\n",
    "        mask = reprojected.difference(cloud_mask)\n",
    "        \n",
    "        if mask.area / reprojected.area < 0.25:\n",
    "            print(f'{product_path}: Area is covered by clouds, skipping it')\n",
    "            return\n",
    "        \n",
    "        miny, minx, maxy, maxx = mask.bounds\n",
    "        \n",
    "        # update the dimensions and save as geotiff, not jp2\n",
    "        out_meta.update({\n",
    "            'width': maxx - minx,\n",
    "            'height': maxy - miny,\n",
    "            'driver': 'GTiff',\n",
    "            'dtype': 'float32'\n",
    "        })    \n",
    "        \n",
    "        out_name = Path(b04_path).name.replace('B04', 'NDVI').replace('.jp2', '.tif')\n",
    "\n",
    "        ndvi_path.mkdir(exist_ok=True, parents=True)\n",
    "        with r.open(ndvi_path / out_name, 'w+', **out_meta) as dst:\n",
    "            # we take only the part out of our source raster that we actually need\n",
    "            # crop=True clips off the borders\n",
    "            b04, transform_b04 = rasterio.mask.mask(b04, shapes=[mask], crop=True, nodata=-999)\n",
    "            b08, _ = rasterio.mask.mask(b08, shapes=[mask], crop=True, nodata=-999) # we ignore the returned transform because it's identical to the previous one\n",
    "\n",
    "            b04 = scihub_normalize_range(b04).astype('f4') # <- f4 = float32\n",
    "            b08 = scihub_normalize_range(b08).astype('f4')\n",
    "\n",
    "            # uncomment the following line to see if your masked shape looks correct\n",
    "            #rplt.show(b04, transform=transform_b04)\n",
    "            #rplt.show(b08, transform=transform_b04)\n",
    "\n",
    "            # we want to be able to ignore divide by zero errors so the formula is nicer to write\n",
    "            np.seterr(divide='ignore', invalid='ignore')\n",
    "            ndvi = (b08 - b04) / (b08 + b04)\n",
    "\n",
    "            # uncomment the following line to see if we calculated the index correctly\n",
    "            # rplt.show(ndvi, transform=transform_b04)\n",
    "\n",
    "            dst.write(ndvi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tqdm.notebook import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "da8f523cda004cf281e368debdddb9c5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(HTML(value=''), FloatProgress(value=0.0, max=40.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "input/tempelhofer_feld/S2B_MSIL2A_20190330T101029_N0211_R022_T33UUU_20190330T144328.zip: Area is covered by clouds, skipping it\n",
      "input/tempelhofer_feld/S2A_MSIL2A_20190626T102031_N0212_R065_T33UUU_20190626T125319.zip: Area is covered by clouds, skipping it\n",
      "input/tempelhofer_feld/S2A_MSIL2A_20190603T101031_N0212_R022_T33UUU_20190603T114652.zip: Problem reading product, skipping it\n",
      "input/tempelhofer_feld/S2A_MSIL2A_20190404T101031_N0211_R022_T32UQD_20190404T174806.zip: Problem reading product, skipping it\n",
      "input/tempelhofer_feld/S2A_MSIL2A_20190216T102111_N0211_R065_T33UUU_20190216T130428.zip: Problem reading product, skipping it\n",
      "input/tempelhofer_feld/S2B_MSIL2A_20190419T101029_N0211_R022_T33UUU_20190419T132322.zip: Problem reading product, skipping it\n",
      "input/tempelhofer_feld/S2A_MSIL2A_20190407T102021_N0211_R065_T33UUU_20190407T134109.zip: Problem reading product, skipping it\n",
      "input/tempelhofer_feld/S2B_MSIL2A_20190512T102029_N0212_R065_T33UUU_20190512T134103.zip: Problem reading product, skipping it\n",
      "input/tempelhofer_feld/S2A_MSIL2A_20190613T101031_N0212_R022_T33UUU_20190614T125329.zip: Problem reading product, skipping it\n",
      "input/tempelhofer_feld/S2A_MSIL2A_20190424T101031_N0211_R022_T32UQD_20190424T162325.zip: Problem reading product, skipping it\n",
      "input/tempelhofer_feld/S2A_MSIL2A_20190822T101031_N0213_R022_T32UQD_20190822T143621.zip: Problem reading product, skipping it\n",
      "input/tempelhofer_feld/S2A_MSIL2A_20190623T101031_N0212_R022_T33UUU_20190623T132509.zip: Problem reading product, skipping it\n",
      "input/tempelhofer_feld/S2B_MSIL2A_20190320T101029_N0211_R022_T33UUU_20190320T195148.zip: Area is covered by clouds, skipping it\n",
      "\n",
      "CPU times: user 191 ms, sys: 172 ms, total: 363 ms\n",
      "Wall time: 12 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "# we parallelize the NDVI calculation using a Python threadpool\n",
    "with Pool() as pool:\n",
    "    for _ in tqdm(pool.imap_unordered(calculate_ndvi, input_files), total=len(input_files)):\n",
    "        # this loop is only here so we can get the progress bar\n",
    "        pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many files could we process?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "28\n"
     ]
    }
   ],
   "source": [
    "! ls -l {ndvi_path} | wc -l"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the NDVI\n",
    "\n",
    "See [02c Multithreading.ipynb](./02c Multithreading.ipynb) for a performance comparison of single vs multi-threaded iteration to average the ndvi."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import ma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_files = list(ndvi_path.glob('*.tif'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "def average(file_path):\n",
    "    with r.open(file_path) as src:\n",
    "        # masked=True makes sure to respect the nodata value and reads the data\n",
    "        # as a numpy MaskedArray, which lets us use `numpy.ma`'s methods\n",
    "        data = src.read(1, masked=True)\n",
    "        return file_path, ma.average(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2f61f096377440319738f5ea049f24b5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(HTML(value=''), FloatProgress(value=0.0, max=27.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "CPU times: user 48.5 ms, sys: 165 ms, total: 213 ms\n",
      "Wall time: 754 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "with Pool() as pool:\n",
    "    averages = [avg for avg in tqdm(pool.map(average, output_files), total=len(output_files))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "paths, vals = zip(*averages)\n",
    "dates = map(lambda p: p.name.split('_')[1], paths)\n",
    "df = pd.DataFrame(\n",
    "    vals, index=pd.DatetimeIndex(dates)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:>"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD2CAYAAAA0/OvUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABAlElEQVR4nO2dd3hc9ZWw3zOakWbU5a7qghsGjG1sAwFCb8kmQCo1ENqSDUtII9n9NrubzWbT2SQLCR1CCBCSQIAEAgQIobqATXHvttybrDqa9vv+uHOl0WgkzUhT7ozP+zx6LM3cOzrHM7rnni7GGBRFUZTDG1euBVAURVFyjxoDRVEURY2BoiiKosZAURRFQY2BoiiKArhzLcBwGDNmjJk0aVKuxVAURckr3n777X3GmLGJnstLYzBp0iSWLl2aazEURVHyChHZMtBzGiZSFEVR1BgoiqIoagwURVEU8jRnkIhgMEhzczN+vz/XogyK1+uloaEBj8eTa1EURVF6KBhj0NzcTEVFBZMmTUJEci1OQowx7N+/n+bmZiZPnpxrcRRFUXoomDCR3+9n9OjRjjUEACLC6NGjHe+9KIpy+FEwxgBwtCGwyQcZFUVxJve8upGX1+zJyGsXlDFQFEUpVCIRw0+eX8sra/Zm5PXVGKSZv/zlL8yYMYOpU6fy/e9/P9fiKIqSJ+xo6eJQZ3DA57cc6KQrGGZWbWVGfr8agzQSDof54he/yLPPPsvKlSt55JFHWLlyZa7FUhQlD7ji3kX8zzOrBnx+9c5WAGbWVmTk96sxSCOLFy9m6tSpTJkyheLiYi6++GKefPLJXIuVNfa0+nn63R10dIdyLYqi5BWdgRAb9nbQ3NI54DGrdrbiEpg+PjPGoGBKS2P59tMrWLmjNa2vOauukv/42FGDHrN9+3YaGxt7fm5oaGDRokVplcPJ3Pf6Zu54ZQMVJW4+Nb+BK06YyJSx5bkWS1Ecz4Y9HQDsbw8MeMyqXW1MHlOG11OUERnUM0gjifZJH07VQy2dASq8bk6fOY6H3trCGT95hc/dt5i/rtxNOKK7thVlINbvbQNgf8cgxmBnK0dmKF8ABeoZDHUHnykaGhrYtm1bz8/Nzc3U1dXlRJZc0NYdYmxFCT+/ZC572o7k0cXb+M2iLVz74FIaanxcccJEPjO/kZqy4lyLqiiOYt3udgAOdgSIRAwuV9+byFZ/kOaDXVyysCljMqhnkEYWLFjAunXr2LRpE4FAgEcffZSPf/zjuRYra7T5Q1R4rTEb4yq83HTmNF77xhncfuk86qt9fO/Z1ZzwvRf5+u/e5YPth3IsraI4h3V7LGMQihha/f0ritbssjyHIzOUPIYC9Qxyhdvt5rbbbuPcc88lHA5z9dVXc9RRufFSckG7P0hFSd+PlKfIxUdn1/LR2bWs3tXKr9/cwuPvbOd3bzczr6maz504ifOPmUCJOzNxUEXJBzbsaafIJYQjhv0dAapL+3rPq6KVRBomyiM+8pGP8JGPfCTXYuSE9u4Q4yu9Az4/c0Il373oGG45byZ/eLuZX7+1hZt/u5z//nMxlyxs4tLjm6it8mVRYkXJPd2hMJv3dzC7oZrl21rY3x7giLhdZKt2tlHl8zBhkL+vkaJhIiVttPlDlJcMfX9R5fNw9cmTefErp/Lg1QuZ01jNbS+v5+QfvMwXHnqbNzfsT5iMV5RCZNO+DiIGjp8yCoADHd39jrGSxxUZLUhRz0BJG+3+EOXe5D9SLpfw4elj+fD0sWw70MlDi7bw2yXbePaDXUwfX84VJ07iE3PrKUvCwChKvrI+mi84YfJo7nxlY7+KokjEsGZXG59d0Jjo9LRRUJ5BPtxN5oOMwyESMbQHehPIqdI4qpR/Of9I3vqXM/nhp2ZT7HbxrT9+wAn/8yL/+dQKNuxtT7PEiuIM1u1uxyVw3KQaoH+vQabHUNgUzC2X1+tl//79jh5jbe8z8HozF/fLFR2BEMbQL4GcKl5PEZ+Z38inj2tg2bYWHnxjM79ZtIUH3tjMKdPGcN0pU/jw9LFDv5Ci5Anr97TTNKqUSq+HCq+bA3GeQabHUNgUjDFoaGigubmZvXszM9EvXdibzgqNNr81gqIihTDRYIgI85pqmNdUw//76Cx+u2QrD721lc8/sIRXbzmdumpNNCuFwfo97UwdZ3Xqjy4r7hcmyvQYCpuCMQYej0e3h+WQ9ug8olRyBskytqKEG8+YxryJNVx69yI27+9QY6AUBKFwhI372jl95jgARpeXsL+9bwI502MobAoqZ6DkjrZoo8xwcwbJUB81ADtadFOcUhhsOdBJMGyYFvUMRpUV9wsTZXoMhY0aAyUt2GGiZEpLh8uEKivXsqOlK2O/Q1GyiV1JNFCYyB5DocZAyRvsMFFlBsJENiXuIsaUl7DzkBoDpTCwjcERtjEotzyDSHSwYzbGUNioMVDSQo9nkEFjAFBf7WW7homUAmHd7jbqq309HvWoshLCMfOJVmdhDIWNGgNlQHYd8nPTI8toPjjwwg2b9iyEiQBqq3waJlIKhvV723u8ArDCRNA7ynplFsZQ2KgxUBLS6g9y1f2LeerdHbyxfv+Qx7f5g4hAWXFmjUFdtY+dLV0F27ynHD5EIob1e9p7ksdghYmgt/Fs9a7Mj6GwUWOg9CMQivCFh95m/Z52RGBHEjH6tu4Q5cXufnPY001dtZeOQJjWLl2tqcDbWw6w7cDQnqsT2d7ShT8Y6WMMRkU9gwMd3T1jKGZOyHyICNQYKHEYY/jGH97j9fX7+cEnZzO2vCSpsIy1yyDzbSt2f8H2JGR6bsUuHn+nOdMiKTlid6ufS+8efIm8k4mvJAIYXVYCwL72AFsPdNIZyPwYCpu0GQMROU9E1ojIehH5ZoLnLxCR90RkuYgsFZGTkz1XyR4/fn4NTyzbztfOmc4nj2ugttrHzkNDJ2xTHVI3XOp6eg2GNgb3v76Jb/7h/aRyHkr+cfvL6+kORXh7y8G8DBuu22NVCk1N6BkEenYYZHoMhU1ajIGIFAG3A+cDs4BLRGRW3GEvAscaY+YAVwP3pHCukgV+s2gLt7+8gUsWNvHF06cCUFflTc4z6A5mtOHMpq7aSqQlU17aFYwQCEf46V/XZVosJctsb+ni0cXbGFVWzJ627qRuWJzG+j3tjK0o6bPIptjt6plPtGpXW1bGUNikyzNYCKw3xmw0xgSAR4ELYg8wxrSbXvNdBphkz1Uyz4urdvOtP37AGTPH8Z0LjupJWNVV+9jR4h/yzqs9yV0GI2VMWQnFRa6kyku7AlZe4fF3mlm7uy3ToilZ5LaXLAP/PxcdDcDybS05lGZ4rNvTztSx5f0eH1Newr72blbtbM3KGAqbdBmDemBbzM/N0cf6ICIXichq4M9Y3kEq514fDS8tdfowunxj2daD3PjwMo6ur+K2S+fiLur9WNRWeekKhjnU1X8vayxt3dnJGbhcwoQkvZWuYJjTZoylrNjNj55bk3HZlOywZX8Hjy1t5tLjmzhj5niK3S6WbT2Ya7FSwhjD+t3tTBvf3xjYIymsSqLs5AsgfcYgUQlJv1tJY8wTxpiZwIXAd1I89y5jzHxjzPyxY3WEcbr4YPshrrxvMeMqS7j3ygWUxpWG1iU5DyhbCWSwQkVJhYkCEWqrfPzjqVN4YeVu3t6SXxcMJTE/e3EdbpfwT6cdQbHbxdF1lXnnGexp66atO9SnkshmVFkxW/Z3su1AdsZQ2KTLGDQDsWt4GoAdAx1sjPk7cISIjEn1XCV9rNnVxhX3LqLC6+E31x7P2IqSfsckm7Bt9w9/sU2q2KGrofAHw/g8RXz+pMmMKS/hB39ZnZeJRqWX9Xva+eOy7XzuxImMizZizW2q4f3thwiGIzmWLnnW7e47hiKWMeXFPdVy2RhDYZMuY7AEmCYik0WkGLgYeCr2ABGZKtFAtIjMA4qB/cmcq6SfDXvbueyeRRS7XTx83fE01JQmPK6uauiEbTAcoSsYzkrOwJLJx65WP6FB/viNMXQFw/iKXZSVuLnpzKks3nSAV9ZqiDGf+elf1+L1FHHDqUf0PDansRp/MNIzxycf2NNm3czUVfUfxW5XFAFZ6zGANBkDY0wIuBF4DlgFPGaMWSEiN4jIDdHDPgl8ICLLsaqHPmssEp6bDrmUxGzd38lldy8CDL+59gQmji4b8Ngx5SV4ioQdg1RrdHRnZxSFTV21j3DEsKet/+Jwm0A4QjhiesJeFy9oomlUKT/8y5qeIWBKX1buaOXlNXtyLcaArN7Vyp/e28nnT5rE6PJeL3ZuUzUAy/IoVGT/zSTa7z0q2mtQ5fNQW5W9rYhp6zMwxjxjjJlujDnCGPPd6GN3GGPuiH7/A2PMUcaYOcaYE40xrw12rpIZtrd0ccndb+EPhXno2uP71DgnIpmEbbq3nA1FMuWl/oDlNdiVGMVuF189Zzord7by9HsahYzHGMNXHlvOFx56u2cCrdP43xfWUlHi5rpTpvR5vL7ax5jykrxKIrd3h4HEN1BjoiMpsjWGwkY7kA8jdrf6uezut2j1B3nomuOTdkFrq3zsHCRGn21jUN/ThTywTF1B64/NF1OW97HZdcycUMFPnl9LIJQ/8eVssGjTAVbvasMfjPD8il25Fqcf7zcf4rkVu7nmlMl96vLBWpE6t6k6r5LIHd0hXAJeT/9LsB0mymaICNQYHDbsa+/msnsWsbetm19dvZCj66uSPreuyjvofCL7TjJbCeTaJJLaPcaguPcj7nIJ3zhvJlsPdPLbpdsGOvWw5FdvbKa61EN9tY8nlm3PtTj9uPWFNVSXerj65MSrbec0VrNxbweHOgcvgXYK7d0hykrcCe/8x1VYnu+sOjUGSppp6Qxw+T2LaD7YyX1XLWBeU01K59dW+9jd6ic8QKzdXnmZrZxBeYmbSq97cGMQ6O8ZAJw2YywLJ43i5y+uozPgzHBIttne0sVzK3Zx8YImPjGvntfX72NPq3M6et/ecpCX1+zl+g9PoXKAGw47b7C8uSV7go2Aju6BmzSnjy/nl5fN44I5dVmVSY1BgdPqD3LFvYvZuK+Dez63gOOnjE75NeqqfQTDhn3tiRO2vZ5BdoyBLdNg5aVdQUsmX1zfhIhwy3kz2NvWzf2vb86kiHnDr9/cAsDlJzRxwZx6Igaeetc5eZVbX1jD6LJirjxx0oDHzG6oRoS8yRt0BEIJk8dgfUbPP6aWEnd2Oo9t1BgUMB3dIT5//xJW72rljsvncfK0McN6nbohdg+3ZmnLWR+ZqgdfctMVTSDHewYA8yeN4qwjx3HHKxto6Qz0e/5wwh8M8+iSrZwzawINNaVMHVfOpNGlvOOQi+pbG/fz+vr9fOG0Iwa8eILlLc4YX5E3eYP27vCg+uQCNQYFSlcgzDW/WsLybS383yVzOWPm+GG/Vm20FnqgYWD2lrOBXPhMUFc9eB4jUQI5lq+fO5P27hC/fGVDRuTLF55cvp2WziBXnTSp57FRZcVDjh/JBsYYbn1+LeMrS7j8hIlDHj+n0Uoi50NjoRUmyu6d/1CoMShA/MEw1/96KYs2HeDWzxzLeUfXjuj16odI2Lb5g7hdQok7ex+numofLZ3BAeP+iRLIscyYUMFFc+t54PXNSY22KESMMdz/+mZmTqjg+Mmjeh6v8nkcsTzotfX7WLz5AF88fWpSw9rmNlXT0hlk837njyzv6A5lfCtgqqgxKDACoQg3PvwOr67bxw8/OZsL5vSb+ZcylT43pcVFA8bo27utXQbZrImuH2Jmkj2xdLCLyJfPmk7EGH7+4uE54npxtJz0qg9N6vPeVfo8OfcMjDH8+Pm11FV5+eyCxqFPAOY0WoUR+ZA3sKuJnIQagwIiFI5w82+X8ddVe/jOhUfz6fnJ/RENhYhQWzXwcLhsDqmzsUNXA3krdjVR/OC9WBpHlXLZ8RN5bGkzG/a2p19Ih/NAtJw0/oah0uuh1Z9bY/DS6j28u62Fm86clnQideq4cspL3HmRN+joDlGmYSIlE4Qjhq/97l2eeX8X//bRI7kiiRhrKgyWsG3zhygvyV6+wJJn8KR2V3DgBHIsN54xlRK3i1ufX5teAR2OXU762QWN+Ir7/h9ZYaJgzmLvkYjh1hfW0jSqlE8e15D0eUUuYXZDFcu2tmROuDTRoQlkJRNEIoZ/ffx9/rh8B18/dwbXxrXrp4O6Kt+A84nau4NZ9wzGV3pxCQPKZOcMhspjjCkv4dpTpvDn93fyXp7UqKeDh96yykkT3TRU+txEDDkbS/Hcil2s2NHKl86chqcotUvU3KZqVu1sxR99/51IIGRt4CvXnIGSTowx/OfTK/jt0m3cdMbUnnWV6aa22svetm66Q/3/yNr8ISqyfJfjKXIxrmLgmUn+YBivx4XLNXQe47pTJjOqrPiwWYDjD4Z5ZHFvOWk8dlWYXTKcTcIRw//+dS1HjC3jwrmp57vmNNYQihg+2H4oA9Klh8GG1OUSNQZ5jDGG7z27mgff3ML1H57Cl8+enrHfZY/a3X2of+NZe5a2nMVTVz2wMegMhAbNF8RS4fXwxdOn8uq6fby+fl86RXQkdjnplR+alPD5Kl/UGOQgifyn93awdnc7N581naIkDHk8cxqrAWevwWzP8pTfZFFjkMf87wtruevvG7nyxIn8y/kzM1rN07PkJkESuc0fymrDmU1dtW/A3oeuQGTIfEEslx3fRH21r+AX4BhjeOCNLcycUMEJU0YlPKYyagyyXVEUCkf46V/XMXNCBR89Znjl0GMrSmio8Tk6b9ARUM9ASSO3v7yen7+0nosXNPIfHzsq42WdtYOMjc7mlrNY6qt9bG/pSnjxtsNEyeL1FHHzWdN4r/kQf/nAeVM708XiTQdYtbO1XzlpLLnyDJ5Ytp1N+zq4+azpSYX3BmJuU42jPYPeMJFWEykj5J5XN/Kj59Zw0dx6vnvRMSP6w0mWuqrEdf3+YNhKhuXgLqe2yksgFGF/R/+REtaWs9T+2D4xr4Fp48r50fNrBt2ils8MVE4ai50zyKZnEAxH+PlL6zi6vpJzjxp+tzxYoaLtLV2OGrYXy2C7DHKJGoM849dvbeG//7yKjx5Ty48+NXtYcdXh4CsuoqbU0y9Gn4shdTaD7WfuDIRSChOBVZr4tXNnsHFvB394pzktMjqJwcpJY6n0We9lNhPIv1vazLYDXXz17Bkj9nKdvvlME8jKiHlsyTa+9ccPOOvIcfz04jm4Uyy7Gym1Vf1j9NlebBNL3SBdyF3BSL+JpclwzqzxzG2q5qd/Xefo8sThMFg5aSx2yC9bYSJ/MMz/vbSOuU3VnDZj7Ihfb1ZtJZ4icWzeQBPIyoh4cvl2vvH4e3x4+lhuv2xeyvXX6SBR9Y49pC7bTWeWPAN7Bv5AGF8KOQMbEWsBzs5D/p7RzoXAUOWksRS5hIoSd9bCRI8u3srOQ36+ds7IvQKw8j+z6qpYvs2ZYynUM1CGzbPv7+Qrj73L8ZNHceflx2V9zrlNoi7ktm7rgpELz6Cm1IPX40poDLqC4ZTDRDYnTBnNqdPHcvvf1ud8LEO6GKqcNJ5KX3ZGUnQFwtz+tw0cP3kUHzoi9V0bAzG3sZr3mg8NuJApl3RGR6VoAllJiZdW7+amR5cxp7Gae69ckHJSNJ3UVvlo9Yd67mygN0yUC5dXRAYsLx1OAjmWr587g5bOIHe9snEkIjqCZMpJ46mMjqTINL9+azN727r5apq8Aps5jdV0BsKs3d2WttdMF+3dITxFkrObuoFQY+BgXl23lxseeocjayu5//MLcu5W1iUoL83FLoNY7PLSeLoC4aTGHg/E0fVVfOzYOu59bRN72pxZlZIsyZSTxlPpdWd8jHV7d4g7XtnIKdPGsHByckYqWXqSyA7MG3Q4cGIpqDFwLG9t3M91Dy5lypgyHrx6Yc4utrHYk0K3xyRse/Yf5yBMZMnUP49hjKErGKZ0hF7UV8+eTjAc4baX1o/odXLNA29spso3eDlpPFVZCBP96o3NHOgI8JUMdM43jSplVFmxI/MG7Q7cZQBqDBzJ21sOcs0DS2ioKeU31x5PdWlxrkUCYjyDmItvrisj6qp97G3vJhDq7QsIhg3hiBl2zsBm0pgyPrugkYcXbWVrjhamRCIm4TyoZNne0sXzK3dz8cLBy0njyfROg0NdQe58ZQNnzhzH3KaatL++iPRsPnMa1pYzNQbKEHyw/RBX3b+YsRUlPHzt8YwuL8m1SD2Mr/QicZNC2/whStwuirO45SyWumofxsDumAYje2LpSMJENjedOQ13kXDrC7kZYvebxVs5+QcvExlmIvSht7ZgjEl5pHmlN7M5g3tf20SrP5TReVpzGqtZt6e9x3t1Ctb4amflC0CNgaNYvauVy+9dRKXXw2+uO4Fxld5ci9QHa1JoSZ+wTFt3bkZR2NT1hK56ZfL3rLwc+R/c+EovV580mSff3cHKHa0jfr1UeXPDPva2dffMs0kFu5z07FnjhywnjafK56EjECaYgU7sgx0B7nttE+cfPYGj66vS/vo2c5uqMQbea3bWBFMnbjkDNQaOYf2edi67exFedxGPXHdCz1pHp2FV7/RNIOeirLRXnv5Lbjp7tpyl5+7rH089gkqvhx89tzotr5cKq3Za1TDD2S3w1PId1rL7D01O+Vy7C7ktA13Id726kY5AZr0CgNkN1YDz1mBqmMgBrNvdxmfueJPvPbuK51bsYm9b/3HMuWDL/g4uu+ctRISHrzueptGp3cVlk7oqHzvjEsi5/GDbjWex5aX2ysuR5gxsqnwevnDaEby8Zi+LNu5Py2smQ2cgxOb9HUBv1VayGGO4/43NKZWTxpKpYXV727p54PXNfGx2HdPHV6T1teOp8nmYOq7ccXkDp1YTOU+iDNLqD9EdjnDfa5u4M2zVjzfU+JjXVMO8pmrmTazhyNrKrHb3Nh/s5NK7FxEIRXj0+hOZMrY8a797ONRWeXlx9W6MMYhIznYZ2Hg9RYwuK+4TJkpnzsDmyhMncf/rm/jhc2v4/Q0nZnxKLMCaXW3YA1nbUvQM7HLS73/imGHJmqlhdXe8soHuUJibz5qW1tcdiDmN1by8ek/P59UJtDvUM3CeRBnkuIk1PPnFk/AHw6zYcYh3trTwztaDLNq0n6fe3QFYaxJnN1Qxr6mGuVEjkanY/a5Dfi67ZxFt/iAPX3cCMyZk9k4pHdRW+/AHIxzsDDKqrJg2f4imUbn1ZGrjxmT05AzSaAx8xUV86czp/OsT7/Piqj2cNWtkkzWTwQ4RAX0a/ZLhV2+mXk4ai73TIJ3lpbsO+XnorS18Yl5D1m565jZV8/u3m2k+2EVjjj+nYHlsHQFnJpAPK2Ng4/UUcdzEURw3sdd93tHSxTtbD7Jsq2Ug7n99M3f+3fIe6qt9zJtoGYa5TTXMqq0ccfXM3rZuLr3nLfa3B/j1NQszmkhLJ/UxMXrbGOQygQxW6MoOp0BvziDd3dqfnt/A3dHx4afPHJfxibGrd/UmrFMJE1nTSXdz7SmTh/1/0BsmGnnOwBjDsx/s4jt/WokxcNMZ2fEKoHfz2TtbDzrCGHSHIoQjprDDRCJyHvAzoAi4xxjz/bjnLwO+Ef2xHfiCMebd6HObgTYgDISMMfPTJVey1FX7qKv28Q+z6wDoDoX5YHsry6IGYunmAzwd4z0cU1/FvIk1zG20wkvjU/AeDnYEuOLeRexs8fPgNQszUmedKezGs52H/BxdX0WbP5jTMBFY792bG3pj+XaYKF0JZBtPkYuvnjOdGx9expPLt/OJeQ1pff14Vu1sZUKll12t/pTCRMMtJ43FTiCPNEy0aV8H//7kB7y6bh9H1lZy26XzspoTmzG+Ap+niOXbWobtJaWTXPflDEZaJBKRIuB24GygGVgiIk8ZY1bGHLYJONUYc1BEzgfuAo6Pef50Y4xjFtCWuIs4bmINx03svVDvPNRleQ5bDvLO1oM88Ppm7oqW3tVX+5gb9RzmNVVzVF1VQu/hUFeQK+5bxMZ9Hdx/1QIWTEpvG36mqY3xDIwxOc8ZgFVR1NYdotUfpNLrwR9If87A5iNH13JM/UZ+8vxaPjq7NmPzZYwxrN7Zxmkzx/H0uzuS9gz8wTCPDrOcNJaqEYaJ/MEwt7+8njtf2UiJ28V/fGwWV5wwMetj191FLo5pqHLMWIqeiaUO7EBOl0QLgfXGmI0AIvIocAHQYwyMMW/EHP8WkNnbqgxQW+Wj9hgfH4nuZ+0OhVm5o5V3oqGld7Yc5E/v7QSgOOo92J7DvKYayr1urrp/MWt2tXHX5+Zz0tQxuVRnWIwpK6G4yMWOQ110BsJETO7vcmJHWVdO8PR4BunMGdi4XMIt583ginsX8/CirXz+pNTLNpOh+WAXbd0hFkyq4el3dySdM3hq+Q4ODrOcNBafpwi3S4blGfx15W7+8+kVNB/s4sI5dfzrR47Mac/M3KZq7n9tM92hcM6Hw7U7dHw1pM8Y1APbYn5upu9dfzzXAM/G/GyA50XEAHcaY+6KP0FErgeuB2hqahqxwOmgxF3E3Gii+RqsP75dh/ws23qwJ//w4FtbuOe1TYAVtugORfjFZfM4fca4XIo+bFwuYUKVl50t/pgtZznOGdjlpS1+Zk6o7DUGGZrwevLUMXzoiNHc9tJ6Pj2/MSPGcNVOK19wdH0VJW5XUn0GIy0njUVEUp5cuu1AJ99+egV/XbWHaePKeeS6EzgxjWOph8vcxmruDEdYuaM15yHZDoeuvIT0GYNEmbSE/fMicjqWMTg55uGTjDE7RGQc8IKIrDbG/L3Pi1kG4i6A+fPnO29IeZQJVV7OP6aW86PeQyAUYeXOVt7ZcpAVO1r5yDETOPPIzFeiZBJ7OFzP+Ooch4nsBj27vNROIHszdBcoItxy3kwuvP117n11E1/KQJnkqp1tiMDMCRVUeN1J5QyWbD7Iqp2tfG+Y5aTxWMPqkvNI3t3WwmfvehOXCP9y/kyuPnlyThYwJcI2AMu3tTjAGNieQeFWEzUDjTE/NwA74g8SkdnAPcD5xpiejJ8xZkf03z0i8gRW2Onv8efnI8VuF3Maq3uqGgqB+mofizYd6Jn5UpHju5wx5SW4XdJTXuoPhvF6XLgyWO0zp7Ga846awN2vbuTyE5rSPkNq1c5WJo0uo7TYTXmJO6mcwQNvbKLK5+HCNCVKK73JbztbvOkA/mCEv33tNCaNKUvL708X4yu91FZ5Wba1hc+flFtZnBwmSpfpXgJME5HJIlIMXAw8FXuAiDQBjwNXGGPWxjxeJiIV9vfAOcAHaZJLyQC11VaFi32hyHUCuSgaurKNQVdg+FvOUuFr506nMxDiF3/bkPbXXr2rlZnRvpNyr3vInMGOaDlpqtNJByOVMJHtjTmhfDMRc5ucMcHUqSsvIU3GwBgTAm4EngNWAY8ZY1aIyA0ickP0sH8HRgO/EJHlIrI0+vh44DUReRdYDPzZGPOXdMilZIbaKh/hiGHjXqu2P9dhIoiu5IyOpBjJystUmDqugk8f18iv39xC88H0jbju6A6x5UAnR9ZWAlblyVBhonSUk8aTkjEIhih2uzLeezFc5jRWs/VAJ/vbczuCpqe01IHVRGkL6hljnjHGTDfGHGGM+W70sTuMMXdEv7/WGFNjjJkT/ZoffXyjMebY6NdR9rmKc7GHw9krBXOdQAaoi/UMgmG8WVoP+qWzpoHAT/+6Lm2vuXl/B8bA1HFWl26Fd/Aw0Uimkw5GpTf5BTf+wMiXCWWS2LxBLrETyE7MGTgjw6PkFXb1jm0MnFAZUVftY9chP+GIoSuLF6a6ah9XnjiRx99pTtu+XTvkYoffykvcg1YTpaucNJ4qn4fWrhDGDF2v0Zml0NxwObquiiKX5LzfoCNg7f/Idr9FMjhPIsXx2F3Ia3e3A84xBqGIYV97d9ZyBjb/dNpUyord/Pi59CzAiR/BPVjOwC4nnTF+5OWk8VT63ATCEfzBoXcadAXDGSvlTQe+4iKOrK3IuWfg1CF1oMZAGQaVXjdlxUXRXa5FjogTx5aXdgXDGek+HoiasmKu//AUnl+5m3fSMDu/K7rIxuexLhplJQPnDOxy0qtOSn7ZfbLYk0uTCRVl0xsbLnMaq3l3W8uwt8alA6eOrwY1BsowEJGeUJET8gXQd0yGP0sJ5FiuPnkyY8qL+cGzq5MKqwxGvGdQUeImEIok3IWc7nLSWFLZaeD0MBHA3MYa2rpDbNjbnjMZ1BgoBUdt1Bg4oZII+o6k6AxkP2RRVuLmn8+YxqJNB3hl7d4RvVa/MFH04mEnH20yUU4aiz3GOpleAytM5IzPwkDMaaoGyGnewAoTOdNoqjFQhkVdlXUn7pT4Z6XXQ0WJmx0tfrqCuQlZXLKwicZRPn74lzUjCkXY+xi8PTkD66IcX1GUiXLSWFIZVmflaZx9OZk8uowqn4dlOcwbdAbC6hkohUVvmMg5H2x7yY0/kN2cgU2x28VXz57Byp2t/On9ncN+nR7PwNPXM4itKMpUOWksld7kx1hbBtg5n4VEuFzCsY3VOd2J3K5hIqXQqI16Bk4yBlbjWVfWms4S8fFj65g5oYKfPL+GYHjoKpxEdAbCFBf1lh8mMgaZKieNpTKFBTedOTLAqTK3sZq1u9tS3hyXLjq6Q45sOAM1Bsow6fEMSpyRQAZLpi37OwlFTM6MgT3iesv+Tn67ZNvQJySgKxDqkwOw8zLt3dYdujGGBzJUThpLTzVRMp5BIOT4aiKw8gYRA+9vP5ST39/RrWEipcCwPQOnJJDBKi+1J6nmsub99BnjWDCphp+9uI6uQP8KoKHojCvTtD0DW7clmw+yMkPlpLEUu134PEVDhomMMTnL06TKnIZqIDdJZGv/sSaQlQKjrtqHp0gYXV6ca1F6sA0U5NYY2COu97Z1c/8bm1I+Pz7MZYfi7GqiTJaTxlPpcw+ZQO4ORYiYzGyWSzc1ZcVMHlPG8m3Zzxt0BsIY48whdaDGQBkmXk8Rv7vhQ1yeoUqW4WCHriAzW85SYcGkUZw5cxy//NsGWjoDKZ3bFVca25szCPaWky7ITDlpPFU+z5CegT9DO6czxZzGapZtbRlxP0iqOHliKagxUEbAnMbqnriyE6h3kDEA+Pp5M2jvDvHLV1IbcR0fJiotLkLEKi21y0mzZYQrvZ4hE8h29ZMT/s+T4ZRpY9jT1s1tL63P6u/tmViqxkBRMsv4Si92CD1bU0sHY+aESi6aU88Dr29mV3S8djJ0xjVwiQjlxW72dQR6ykmztTfA2nY2uGfQYwwc8H+eDBfNrecTc+v5yQtreWJZc9Z+b+/EUjUGipJRit0uxkY3jpU65C71y2dPJ2IMP3sx+RHXXYFQP/nLvW6eftcqJ73yQ5PSLOXAVKYUJnLmRS4eEeH7n5zNiVNGc8vv3+OtjfuHPikNtDt45SWoMVAKDDtv4JS71MZRpVx2/EQeW7qNjUnOxIkPE4EVWmjzh5gxvoITp2RvyXyl1z1kaWm+hYnAunG44/LjmDi6jOsfXMr6PekZPz4YHRomUpTsYecNnHRh+uLpUylxu/jJC2uHPpjoDuc4Y2CHFjJdThpPlc9DW3do0PEaXcH8ChPZVJV6uP+qBRS7XVx1/xL2tmV2C1pHQBPIipI17PJSJ5U5jq0o4dqTJ/Pn93byfvPQzU6dgXC/MFGlz5O1ctL432sMg67d7B257Zz/82RpHFXKvVcuYF97N9c+uHRYfSHJoglkRckiU8aW43ZJzygFp3Ddh6fg9bj4wzuDJywHauD6+jkzuPOK47J+912ZxBjr+Cmr+caxjdX8/OK5vNfcws2/XUY4Q/sOtLRUUbLIp45r4M83ndIzcdMpVHg91JQWD7q+EsAfjGAM/cZBH9NQxQlZzBXY2KXDgyWRu/KszyAR5xw1gW99dBbPrdjN955ZlZHf0d7ddwCh03CmiVKUYVLsdjFjQkWuxUiI11PUU3kzEJ3RkItTLqyVPusSMVh5qR1acUI570i4+uTJbD3QyT2vbaJxVGnaq7Y6opsBXQ7YDJgINQaKkiVK3K4h9wn3JGMdcvdYlcTk0q48rCYaiG/9wyyaD3bx7adXUF/t46xZ49P22k7ecgYaJlKUrOH1FCVcXRlLl8MauJKZXNoZDOMpEjxF+X85KXIJP79kDkfXV/HPjyxLKuGfLNaWMzUGinLY4/W4kggTOSv+XpnEtrOuPNh/nAqlxW7uuXI+o8qKufpXS2g+2JmW11XPQFEUwM4ZDB4mctpoh4oSNyJDJJBzsHM604yr8PLA5xfgD4a5+oElSW17Gwprl4Fz/5/UGChKlvC6h04gdwXtBLIz7iBdLqGiZPAu5M48WHk5HKaNr+DOy49j494OvvDQ2wRCw9tcZ9PeHaLMwf9PagwUJUt4PS78Q+YMrAuOU8JEYHXqtvoHTyAXUpgolg9NHcP3PzmbNzbs51+feH9EY687As4OEzlXMkUpMLyeIrqHDBM5r5u30jv4sLquYKjgwkSxfOq4BrYd6ORnL66jaVQpN505bViv4/ScgXMlU5QCI5k+AyfO+bF2GgyeM3DyRS4d3HzWNLYd6OTWF9bSOMrHRXMbUn4Nq5rIOe9rPBomUpQsUeJx4R8i7uy0aiIYettZZyDsqFlQmWCkY69D4Qj+YMTRRlONgaJkCa+7iEAoMugEUNsYeN3OubgOtQc50SylQmQkY687ox7fYdFnICLnicgaEVkvIt9M8PxlIvJe9OsNETk22XMVpRCw7567B/EO/EErGeukkQVVvsFXX3Yl2L9QqAx37LXTh9RBmoyBiBQBtwPnA7OAS0RkVtxhm4BTjTGzge8Ad6VwrqLkPV6P9ec2WN6gM+C8ZGyl10NXMDxgaWXXYRAmiqVxVCn3pDj2+rAxBsBCYL0xZqMxJgA8ClwQe4Ax5g1jzMHoj28BDcmeqyiFQEk09DNYeWmnA8s0h+pCPlzCRLHMaazmZymMvbYnlh4OCeR6YFvMz83RxwbiGuDZVM4VketFZKmILN27d+8IxVWU7NPrGQwcJnJiyKVqkJ0GgVCEUMQ4zoBlg3Njxl7/zxBjr3s8Awc3naVLskQBzoSmUkROxzIGJ6dyrjHmLqKhpfnz52dm+4SiZBA7lDJ4mMh5xsAeY52ooqh3sJ5zL3KZxB57fe9rm2is8XHVSZMTHteeB2GidEnWDDTG/NwA7Ig/SERmA/cA5xtj9qdyrqLkO8nkDJw456dncmmCLuRCWGwzUuyx1//1p5XU15RydoKx1x0OX3kJ6QsTLQGmichkESkGLgaeij1ARJqAx4ErjDFrUzlXUQoBu1x00DBR0Hk5AztMlMgzcGLHdLaJHXt90wBjrw+bBLIxJgTcCDwHrAIeM8asEJEbROSG6GH/DowGfiEiy0Vk6WDnpkMuRXESJZ5kEsghxw19G2wPshM7pnPBUGOvexPIznpvY0lbn4Ex5hljzHRjzBHGmO9GH7vDGHNH9PtrjTE1xpg50a/5g52rKIWGHSbqzrMwUdUg1USFtOVspIyr8HL/AGOvO7pDuKT3M+BEnCuZohQYvQnkgcNEnQ4s0yxxuygucg0QJtKcQSzTBxh73R4dUifinGbCeNQYKEqWSLaayGmegYhYIykSdCHbYaLDqelsKBKNve5w+MpL0KmlipI1vO7Bq4nCEUMgFKHU47w/y0pf4smlXeoZJCR+7LXTdxmAGgNFyRpDzSbqTcY6z2Gv9HoS5wx6Skv1UhJP7Njr6lIPE0eX5VqkQXHep05RCpShcgY9ZZoOvLBWDeAZdGoCeUDssdcnTBlFS2fQ0aMoQI2BomSNIpfgKZIBS0t7Qi4OvLBWDrDToKvHgDlPZidQ7HZx5+XzmTmhgiljynMtzqA47xZEUQoYr3vgbWdOrsyp9LoH7EC2jZySmKpSD3++6RQcNJU8IeoZKEoWKfEUDRgmcnIDlx0mil8I3xkIU+opcnTJpBMoconj/4/UGChKFvF6XAM2nTm5gavS5yEUMT3ei40/6LxSWGV4qDFQlCzi9RQNmDPoDRM5L3o7UBeyE/silOGhxkBRsojX40qimsh5F1d7cml8ErnLgct4lOGhxkBRsshgCWQnN3DZOw3iu5C7NExUMKgxUJQs4vXkZzXRQNvOnLiMRxkeagwUJYsMFiZycjXR4GEi5+U4lNRRY6AoWaRkkARyVyCMS6C4yHl/lpUDJJA1TFQ4OO9TpygFjNddRPeACeQwpcXOHHNc6U28B7kr2meg5D9qDBQli1hhogE8g2DIsXfZ7iIXZcVF/RLInQHnyqykhhoDRckiQyWQnZyMrfL1n1yqYaLCQY2BomQRr8c18Ahrh9fsxw+rC4YjBMPG0TIryaPGQFGySIm7iFDEEAr3NwhOv8uu9PYdY927y8C5MivJo8ZAUbKIvRDdn8A7cHqYqNLn6TO51B9wbimskjpqDBQliwy2B7nT4TX71h7kXs9AF9sUFmoMFCWLeN0DG4OuQMjRnkH8tjMnd0wrqaPGQFGySIkdJkrQa9AVdHiYyOuhrTtEOGLtNLBzBl71DAoCNQaKkkWGChM5+cJqdyG3RctLuxw8cltJHTUGipJF7It9d4KRFF0OTyD3DquzkshaTVRYqDFQlCzidScOEwVCEUIR4+gLa/xICnv/gpO9GSV51BgoShYZKEzUs/LSwSGX+GF1Tt6/oKSOGgNFySK9xqCvZ5APIRc7TGR7Bj0jt9UzKAjUGChKFrEv9h3d/Qe+xT7vRCrjFtx0atNZQaHGQFGyyPhKLyKw41BXn8ftC6uT4+9VcWEif9Dav1Di1stIIZC2d1FEzhORNSKyXkS+meD5mSLypoh0i8jX4p7bLCLvi8hyEVmaLpkUxWkUu11MqPTSfLCvMciHMFFZcREuiU0gW4P1nLh/QUmdtGSrRKQIuB04G2gGlojIU8aYlTGHHQBuAi4c4GVON8bsS4c8iuJkGmp8NB/s7PNYPnTziog1nyhaWtoZCDs64a2kRro8g4XAemPMRmNMAHgUuCD2AGPMHmPMEiCY6AUU5XChoaa0v2cQzRk4eTYRWKEi2zPwB8P4ijVEVCik652sB7bF/NwcfSxZDPC8iLwtItcnOkBErheRpSKydO/evSMQVVFyS0ONj52H/H3GWOdDmAiiY6z9vX0GpQ43XkrypMsYJAoamhTOP8kYMw84H/iiiHy434sZc5cxZr4xZv7YsWOHK6ei5Jz6ah/hiGFXq7/nsXypzIkdVtcVjOB1uLxK8qTLGDQDjTE/NwA7kj3ZGLMj+u8e4AmssJOiFCQNNaUAfUJFXXliDCp97t4+g0CIUgdXPympkS5jsASYJiKTRaQYuBh4KpkTRaRMRCrs74FzgA/SJJeiOI6GGh/Q1xj0JJAdfnG1wkS9s4mcHtZSkictAT9jTEhEbgSeA4qA+4wxK0Tkhujzd4jIBGApUAlERORmYBYwBngiWp7mBh42xvwlHXIpihOprbZ6DWIrijoDYYqLXLiLnJ2QjQ0TdQbCGiYqINKW/THGPAM8E/fYHTHf78IKH8XTChybLjkUxemUuIsYX9G318Dv8P3HNpU+D92hCP5g2Jqy6nBPRkkeZ9+GKEqBEt9r0OnwLWc29uTSVn+QrjwxYEpyqDFQlBxgGYO+OYN8GPgWO5/IajpzvsxKcqgxUJQc0FBT2qfXoCtPLqy2MTjYGSQQimifQQGhxkBRckBDTd9eg06HbzmzsYfV7Y7KrR3IhYO+k4qSA+J7Daz4u/Pvsiu9ljHYdcg2Bs6XWUkONQaKkgPiew3ypTKn0mdd/Hs8gzyQWUkONQaKkgPiew06g6H8yBnYnkFrN+D8WUpK8qgxUJQcEN9rkC8JZK+niBK3i909YSLny6wkhxoDRckRsb0GnXkSJgIribxLw0QFhxoDRckRdq+BMSav5vxUxhiDfJFZGRo1BoqSI+xeg45AGGPypzKn0usmELL6I9QzKBzUGChKjrB7DTbt7QDA58mPP0e71wA0Z1BI5MenT1EKELvXYO3uNgBK88UziDUG6hkUDGoMFCVH2L0Ga/dYxiBf7rLt8lLIHwOmDI0aA0XJEXavwbrd7UD+JGNjw0TePAltKUOj76Si5Ai712BdvnkG0S5kn6eI6FIqpQBQY6AoOaShxse2A1bjWb6EXGzPIF88GSU51BgoSg6x8waQP8lYO2fgzRN5leRQY6AoOcSuKIL8udOuVM+gIFFjoCg5pI9nkCcXVw0TFSZqDBQlh+SlZ6BhooJEjYGi5JBYz8Drzo+Lq3oGhUl+lC8oSoFi9xqUuF24XPlRplnujZaWqjEoKNQYKEoOsXsNAuFIrkVJmiKXUFHixufRy0choe+mouSYhhofO6PLYvKFy0+cyJzG6lyLoaQRNQaKkmMWTh7F2uhIinzhG+fNzLUISppRY6AoOeYWvbAqDkCriRRFURQ1BoqiKIoaA0VRFIU0GgMROU9E1ojIehH5ZoLnZ4rImyLSLSJfS+VcRVEUJbOkxRiISBFwO3A+MAu4RERmxR12ALgJ+PEwzlUURVEySLo8g4XAemPMRmNMAHgUuCD2AGPMHmPMEiCY6rmKoihKZkmXMagHtsX83Bx9LNPnKoqiKGkgXX0GiYaqmHSeKyLXA9dHf2wXkTXR78cA+5L8XflCIepkU6i6FapeNoWsX6HqlkiviQMdnC5j0Aw0xvzcAOxI57nGmLuAu+IfF5Glxpj5yYvqfApRJ5tC1a1Q9bIpZP0KVbdU9UpXmGgJME1EJotIMXAx8FQWzlUURVHSQFo8A2NMSERuBJ4DioD7jDErROSG6PN3iMgEYClQCURE5GZgljGmNdG56ZBLURRFSY60zSYyxjwDPBP32B0x3+/CCgEldW4K9AsdFQCFqJNNoepWqHrZFLJ+hapbSnqJMcnmeRVFUZRCRcdRKIqiKGoMFEVRFDUGiqIAIpIfC5iVjJEXxqAQP6giMirm+4LRT0ROE5GxuZYjE4jIV0XknOj3BfOeRamwvyk03QpNn1jSqZujjYGIXCAivwKOzbUs6SI6ofXvwE9F5CcApgCy+DF6XQZ051qedCIi54jIc8A3gM9BYbxnACJytoi8BvxYRG6BgtKt4K4fNpnQzXFrL0VEjDFGRE4HvoM12O5EEdlijDmYY/GGRdR6u4BrgKuB7wHLgAdF5HxjzLO5lG+4RPUS4LPAncA1xpjf5Vaq9BDVzQP8O3Aq1ntWDCwQEQ8QyveLpog0AP8JfB/4G/CoiIw2xnzD/jvMpXwjoZCuHzaZvjY6yjOI+wBuAs4Fvg4cD8zOmWAjwNbJGBMGXgNONsY8CfiBPcAKEXHZx+ZQ1JSI0SuCNT7kQWB99LlPi0hD9KKZV3pBH90CwJPGmFOivTAHgYuNMcF8vVDGvRczgfeNMU8bY9qwRsl/WUSmRS86efW+xbEJOIc8v37YZOPa6BhjEO1CflxEviwiE4wxm40xO40xLwG7gVNFJK+mmcbpVGuMWRnt1p4H/BGYhBV6uNU+JTeSpkaMXl8RkTFYRu494Jcishr4DPB/wC/sU3IjaeokeM+WRB/3GGNeATaKyPm5lXJ4xOlWCawFThaRE6OHjANWAP+WKxmHi4j8k4h8Mvq9ANuMMbvy+fphk61royOMgYhcBFwJ/BzLyv2biMyJOeQ3wHQsKxh7nmMvMgl0+n8xOtl3mAuBW4CrRGR+9C7b0cTpdQzwbWAq8CfgZeASY8ynsUJiF4rIcfmgFwz4ntkx2VA06b8FCOdIxGGTQLcfYOV2/hf4RxF5HetO+hPAHBGZlA/ej4hUiMgdWOG8X4mIOyp3rGeTd9cPm2xeGx1hDLAU+aUx5mWsGOYmrK1oABhj3sMaaHe0iJwhIt+IPu7kD2sinb4EYIzZZIzZGv2+A3gMa2ZTPhCv12bg68aYHcC3jTHLAIwxB7C8n/LciDksBnvPTFQnH3A6gB3eyxMS6fZtY8y9wHXAl40xlwJbgcVAa64ETYVoeOsVY8wErBuS26NP9YRV8vT6YZO1a2NWP8zx1irm543ApQDGmC3An4EyEfl4zOGPANcCv8Wa0+0Iy56iTqVxOiEi/wYcBazMvLTJk4JeTwMVIvJxY4w/5vhvYem1OjsSJ88IP4cPAQtFxOtEjycF3Z4CakTkomgOZHH0uO8AZUBblkROmkF0s6cc34y1NneaMSYsIu6YYxx5/bBxwrUx23c2faqXYqzX74FOEbHXXe7Eqm6YJRblwM+A94HZxpivx52fS1LWCUBEzherpG868KnoID8nMVy9ThGRl7H0+qQxZnd2xE2JYX0Oo4/5sFazOjVUlKpuMwBEZJqIPAkcjeUlxK+ndQIJdTPGdIiIK/o39AvgnujjoWgivAwrzOLE64eNJ/aHXFwbs2IMROQEEfkN8O3oh64o+rj95h4EngC+ICJijDmEFV7wRpXyA18yxnzUGLMzGzIPxQh08kWfXwXcYIz5nFN0grTotRn4ojHmCifpBSPSrSTmj+tJY8zdTrtYjuRvLPr8Lqz37eNOM+CD6FYUH6ozxnwTmCwiJ4rIeBFZEA3F3uSk64dNVM7fAT8SkVm5vDZm3BiIyNFYlSV/wiqlvJ7exp1Q9DAf1j6DncBdIlIHzMWqo7Ut/J5My5osI9QpED1uszHmgyyLPihp0mubMcZRIS8YsW728xirRNhRpOlvrM0Y05xl0YdkCN3CxphI9O64Kua0HwCvA68CpdFjHXP9sBGRccBtWOP792Plp66G3Fwbs+EZnACsNsY8AtwNdAKXicgUABH5DpblGw98FatU6mGgBasZxokUok5QuHqB6lbIuv0eK7yFWGW//4xVrn2UscqBncqxwFpjzP3AT4DHgQtEZCaAiPw32XzfjDFp/cLq1jw+5udjsWJcU6M//wfWm/dtrETVw8ARca9Rmm65VKfDRy/V7fDVDSt31ZhrPQbQ7ULgX4GPRn8eC6yz5QdGRfX7AZZHk9X3LW3jKESkAvgVcBrwRxFZZ6xSvA1YpWr3icgBrCTQQ8B8oMtY5WxEE0ARAGNMZ7rkGgmFqBMUrl6guh3GuhUZK2zkxBDlWCyvZhTWBf5+EfknY8zvReQPWJ7MzVh3/C9i9RV4s/2+pTNMFABeAi7HGk/waQBjTLsx5hbgRuB+Y8w/YI0tOMpWMFZZh1GIOkHh6gWq2+Gqm+NyOTEcAbxujPmwsVYBfxX4SvS5R4CZInJWVJ/9WGGhbsju+zYiz0BEPofVkfmuMaZFRO4BIli1rieLyHRjzFroaY54L3rqGcBb0ey4cdKHtBB1gsLVC1Q31c2xutkNfG9jNYsRrRZaiTX2A6yS0EexphhfCJyJNb7FA5BN3VLegSwiAkzAcnciWG5cGVZ5077oMdOwXB2/Mea/Y849DitREgauN8ZsSIcSI6UQdYLC1QtUN9Ut/3SzQ1kicjnwcWPMZ2LOvQWrN2cmcJ0xZlXWFTCpJUCKov9OBx6Kfu/GKv36Q9yxF2E1gEwFfNHHRgOnpvI7M/1ViDoVsl6qm+qWaz2Godvjccc8CHwm+v2EmNcozqUOSYWJxGqA+C+gSESewZqjEwarzlVEbgJ2iMipJlrKZYx5QkSOBP4ClIvIGcZK7jii1KsQdYLC1QtUN1S3gtANaAc2ich/AZ8QkfOMMc3GGpmeM4ZMIIvIqVgxrxqsxI29VOF0EVkIPa3P/4U1SMk+79PA/8OaZDnbOCjLX4g6QeHqBaqb6lYYukVzBldjlcdWAqcbpzT7JeH+nAJcEfPzL4AvAFcBb0cfc2HFyh4DJsecd0ou3Z7DSadC1kt1U92c+DUM3SZiVRb9FJiXa/njv5IpLX0beCxq0cBq824yxjyA5Rr9s7Ey3g1A2BizCcAY86ox5tUkXj8XFKJOULh6geqmujmPVHSLGGO2GGM2GGNuNsa8kyOZB2RIY2CM6TTGdJveOt6zgb3R7z8PHCkif8Kql3WcgokoRJ2gcPUC1U11cx4p6vY2OGtsdjxJ9xlErZ/Baoiw54e3YbVXHw1sMsZsT7uEGaQQdYLC1QtUN9XNeaSim4nGjpxIKh3IEaxGiH3A7KjF+xaW+/Nanr6RhagTFK5eoLrlK6qb00klwYA1QTCCtQD9mlwlOtL5VYg6FbJeqlv+fqluzv5KqQNZRBqAK4BbjTHdwzM/zqIQdYLC1QtUt3xFdXM2KY+jUBRFUQqPbO9AVhRFURyIGgNFURRFjYGiKIqixkBRFEVBjYGiKIqCGgNFURQFNQaKoigK8P8B9I7yjpFMgNoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}