tweets.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "import time\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import folium\n",
    "import gc \n",
    "from IPython.display import HTML\n",
    "import warnings\n",
    "import geopandas as gpd\n",
    "from shapely.geometry import MultiPoint, Polygon, Point, LineString\n",
    "from sklearn.cluster import MiniBatchKMeans, DBSCAN\n",
    "from sklearn import metrics\n",
    "import sqlalchemy\n",
    "import psycopg2\n",
    "import os\n",
    "import subprocess\n",
    "from scipy.spatial import ConvexHull\n",
    "# from scipy.spatial import Delaunay\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "%matplotlib inline  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup database connection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "USER = 'harry'\n",
    "PASSWORD = 'harry'\n",
    "enginestring = 'postgresql://{}:{}@localhost:5433/tweets'.format(USER, PASSWORD)\n",
    "engine = sqlalchemy.create_engine(enginestring)\n",
    "\n",
    "conn = psycopg2.connect(\"dbname=tweets user=harry host=localhost port=5433 password='harry'\")\n",
    "cur = conn.cursor()\n",
    "conn.commit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"tweets_1M.json\",'r') as f:\n",
    "    tweets=json.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1000000 rows of data\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>user_id</th>\n",
       "      <th>text</th>\n",
       "      <th>timeStamp</th>\n",
       "      <th>lat</th>\n",
       "      <th>lng</th>\n",
       "      <th>id</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>224874450</td>\n",
       "      <td>@Tanner_Cortez hey checkout the website: http:...</td>\n",
       "      <td>Wed Sep 11 04:38:08 +0000 2013</td>\n",
       "      <td>37.446100</td>\n",
       "      <td>-121.883557</td>\n",
       "      <td>377652254096228352</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>312179473</td>\n",
       "      <td>@Lennayy7 😪</td>\n",
       "      <td>Wed Sep 11 04:38:08 +0000 2013</td>\n",
       "      <td>34.087406</td>\n",
       "      <td>-117.462604</td>\n",
       "      <td>377652255346159616</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>54351774</td>\n",
       "      <td>i laugh a lot with that line</td>\n",
       "      <td>Wed Sep 11 04:38:10 +0000 2013</td>\n",
       "      <td>37.356131</td>\n",
       "      <td>-121.842867</td>\n",
       "      <td>377652262325456897</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>343219606</td>\n",
       "      <td>sons of anarchy is back on woop woop</td>\n",
       "      <td>Wed Sep 11 04:38:11 +0000 2013</td>\n",
       "      <td>37.364664</td>\n",
       "      <td>-122.009629</td>\n",
       "      <td>377652264682655744</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1569395935</td>\n",
       "      <td>Drinking a Fresh Squeezed IPA by @deschutesbee...</td>\n",
       "      <td>Wed Sep 11 04:38:12 +0000 2013</td>\n",
       "      <td>37.382600</td>\n",
       "      <td>-121.995000</td>\n",
       "      <td>377652271116722176</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      user_id                                               text  \\\n",
       "0   224874450  @Tanner_Cortez hey checkout the website: http:...   \n",
       "1   312179473                                        @Lennayy7 😪   \n",
       "2    54351774                       i laugh a lot with that line   \n",
       "3   343219606               sons of anarchy is back on woop woop   \n",
       "4  1569395935  Drinking a Fresh Squeezed IPA by @deschutesbee...   \n",
       "\n",
       "                        timeStamp        lat         lng                  id  \n",
       "0  Wed Sep 11 04:38:08 +0000 2013  37.446100 -121.883557  377652254096228352  \n",
       "1  Wed Sep 11 04:38:08 +0000 2013  34.087406 -117.462604  377652255346159616  \n",
       "2  Wed Sep 11 04:38:10 +0000 2013  37.356131 -121.842867  377652262325456897  \n",
       "3  Wed Sep 11 04:38:11 +0000 2013  37.364664 -122.009629  377652264682655744  \n",
       "4  Wed Sep 11 04:38:12 +0000 2013  37.382600 -121.995000  377652271116722176  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.DataFrame(tweets)\n",
    "print (len(df),\"rows of data\")\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convert to state plane coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['lat_conv'] = (df.lat - 32)*89700 # convert to meters with origin at 32d N (1d lat = 89.7km)\n",
    "df['lng_conv'] = (df.lng + 126)*112700 # convert to meters with origin at 126d W (1d lng = 112.7km)\n",
    "df = df.assign(crds=df[['lat_conv', 'lng_conv']].values.tolist())\n",
    "# # define the number of kilometers in one radian\n",
    "# kms_per_radian = 6371.0088\n",
    "# # define epsilon as 1.5 kilometers, converted to radians for use by haversine\n",
    "# EPSILON = 1.5 / kms_per_radian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run Mini-Batch KMeans to divide data into spatial clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 10 # set the number of clusters to split the data into\n",
    "def runMBKmeans(df):\n",
    "    # RUNNING MINIBATCH KMEANS ON ENTIRE SET \n",
    "#     print ('Beginning to run MiniBatch Kmeans.')\n",
    "    ## initialize with MBK-means++, a good way to speed up convergence more than K-Means\n",
    "    mbk_means = MiniBatchKMeans(init='k-means++', n_clusters=n, n_init=10, init_size = 3*n, verbose = 0) # set minibatch kmeans parameters\n",
    "    mbk_means.fit(list(df.crds.values)) # run mini batch k means on the lat and lng coordinate data\n",
    "    df['mbk_means_labels'] = mbk_means.labels_ # this is an array that indicates which cluster each tweet is in\n",
    "    df_output = df\n",
    "    mbk_means_cluster_centers = mbk_means.cluster_centers_ # this is an array that indicates which cluster each tweet is in\n",
    "    df0 = pd.DataFrame()\n",
    "    df0= pd.DataFrame(mbk_means_cluster_centers,columns = ['lat','lng'])\n",
    "    df0 = pd.DataFrame(mbk_means_cluster_centers,columns = ['lat_centroid','lng_centroid'])\n",
    "    df0['unique_tweeters'] = df.groupby('mbk_means_labels')['user_id'].nunique()\n",
    "    df0['unique_tweets'] = df.groupby('mbk_means_labels')['mbk_means_labels'].count()\n",
    "#     print ('Finished running Kmeans.')\n",
    "    return df_output, df0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "df,df0 = runMBKmeans(df)\n",
    "df0['labels'] = np.arange(n)\n",
    "df['db_labels'] = \"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running DBSCAN on each cluster generated by MiniBatch Kmeans"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def runDBSCAN(df,EPSILON):\n",
    "    coords = list(df.crds.values)\n",
    "    db = DBSCAN(eps=EPSILON).fit(coords)\n",
    "    db_labels = db.labels_ # create DBSCAN cluster labels generated for the cluster generated by MBKmeans\n",
    "    db_labels_unique = np.unique(db_labels)\n",
    "    k_DB = len(db_labels_unique)-1   #number of clusters detected with DBSCAN, subtract one because one cluster is noise (-1)    \n",
    "    DB_noise = list(db_labels).count(-1) # these are tweets that DBScan determined to be noise\n",
    "    df.db_labels = list(db.labels_)\n",
    "#     print ('Number of sub-clusters detected is: ', k_DB)\n",
    "#     print ('Noise tweet count is: ', DB_noise)\n",
    "#     print ('Important (non-noise) tweet count is : ', len(df)-DB_noise)\n",
    "    collected = gc.collect() \n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cluster # 0 Run Time:  13.0 seconds\n",
      "Cluster # 1 Run Time:  5.0 seconds\n",
      "Cluster # 2 Run Time:  2.0 seconds\n",
      "Cluster # 3 Run Time:  1.0 seconds\n",
      "Cluster # 4 Run Time:  6.0 seconds\n",
      "Cluster # 5 Run Time:  2.0 seconds\n",
      "Cluster # 6 Run Time:  29.0 seconds\n",
      "Cluster # 7 Run Time:  2.0 seconds\n",
      "Cluster # 8 Run Time:  2.0 seconds\n",
      "Cluster # 9 Run Time:  7.0 seconds\n"
     ]
    }
   ],
   "source": [
    "clusters = np.arange(n)\n",
    "EPSILON = 400 #meters\n",
    "for cluster in clusters:\n",
    "#     print ('________________________________________')\n",
    "    start_time = time.time()\n",
    "    df_cluster = df[df['mbk_means_labels']==cluster]\n",
    "    if len(df_cluster)>1:\n",
    "        df_cluster1 = runDBSCAN(df_cluster,EPSILON)\n",
    "        df[df['mbk_means_labels']==cluster] = df_cluster1\n",
    "        elapsed_time = time.time() - start_time\n",
    "    else:\n",
    "        next\n",
    "    print('Cluster #', cluster, 'Run Time: ', round(elapsed_time,0), 'seconds')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Removing points determined to be noise in DBSCAN (i.e. cluster labels of '-1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Rows of data reduced from 1000000 to 985158.\n"
     ]
    }
   ],
   "source": [
    "df = df[df.db_labels!= -1]\n",
    "print ('Rows of data reduced from 1000000 to' , str(len(df)) + \".\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a unique cluster label that combines both kmeans and dbscan cluster labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3545 clusters were found in the data.\n"
     ]
    }
   ],
   "source": [
    "df['labels'] = (df.mbk_means_labels).astype(str) + '.' + (df.db_labels).astype(str)\n",
    "num_of_clusters = len(df.groupby('labels').labels.nunique())\n",
    "print (num_of_clusters , \"clusters were found in the data.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Write data to  psql / postgis database"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.to_sql('tweetpoints', engine, if_exists='replace')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "cur.execute(\"SELECT AddGeometryColumn ('tweetpoints', 'location', 4326, 'POINT', 2);\")\n",
    "cur.execute(\"UPDATE tweetpoints SET location = ST_SetSRID (ST_Point(lng, lat), 4326);\")\n",
    "conn.commit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a new dataframe to hold cluster information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1 = pd.DataFrame()\n",
    "df1['label'] = df.groupby('labels').labels.nunique().index.values\n",
    "df1['tweetcount'] =  df.groupby('labels').user_id.count().values\n",
    "df1['usercount'] =  df.groupby('labels').user_id.nunique().values\n",
    "df1 = df1.sort_values('tweetcount', ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>label</th>\n",
       "      <th>tweetcount</th>\n",
       "      <th>usercount</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2591</th>\n",
       "      <td>6.0</td>\n",
       "      <td>241859</td>\n",
       "      <td>27035</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.0</td>\n",
       "      <td>164483</td>\n",
       "      <td>9278</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3343</th>\n",
       "      <td>9.0</td>\n",
       "      <td>85945</td>\n",
       "      <td>8091</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1950</th>\n",
       "      <td>4.1</td>\n",
       "      <td>73890</td>\n",
       "      <td>3655</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1949</th>\n",
       "      <td>4.0</td>\n",
       "      <td>41014</td>\n",
       "      <td>6676</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     label  tweetcount  usercount\n",
       "2591   6.0      241859      27035\n",
       "0      0.0      164483       9278\n",
       "3343   9.0       85945       8091\n",
       "1950   4.1       73890       3655\n",
       "1949   4.0       41014       6676"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df1.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Filter out clusters that contain less than a thousand tweets in them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "186 clusters contain at least 250 tweets.\n"
     ]
    }
   ],
   "source": [
    "threshold = 250\n",
    "df2 = df1[df1.tweetcount >= threshold]\n",
    "df2 = df2[df2.usercount > 1] # remove clusters that consist of only one user account\n",
    "print (len(df2), \"clusters contain at least\", threshold,\"tweets.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Investigate distribution of the number of tweets in each cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA70AAAGMCAYAAAAWZMWNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de7y1c53/8dc7yqgochdDupGazqo79KuEIknTScV0oDRSmabTL53TNJ3MSDVKUQaNHDqbMkWOv5rITU45EyXi5hYihM/vj+vaWVZ73/dee6+999qX1/PxuB57re/1Xdf1Weu+rO29r+/1vVJVSJIkSZLURfeb6wIkSZIkSZophl5JkiRJUmcZeiVJkiRJnWXolSRJkiR1lqFXkiRJktRZhl5JkiRJUmcZeiVJGkeSPZNUkoVzXcuyJFnY1rnnXO5zLuqYy/1KkuYPQ68kaUYl2bwNJWPLXUluSHJukoOTbJMkc13nZLTvZc8kD52h7fd+TncnuTnJZUm+m+QNSVYe8v52TvKOYW5zJrTBds8kG811LZKk+SdVNdc1SJI6LMnmwAnAYcDRQIBVgMcCLwXWBX4CvLKq/jBHZf6VJCsCKwK3V/vLsj2b+FFgvaq6fAb2WcCZwN5t0wNpPp+tgWcAlwKvqKqzel4TYCXgzqq6c8D9nQgsrKqFA77ur/bZnhH/NfCxqtpzkO1NYn+b0xxDb6iqg5ZXiyRJvVac6wIkSfcZZ1TVf/U2JHkXsBfwLppQ/MK5KGw8bYCaixD1u/7PCfhQklcChwL/k+QJVXUDQBvIb5uNwpKsUlU3z+Y+l2eUapEkjSaHN0uS5kxV3VVV7wZ+CmyT5Nm965M8JMlnklyS5PYkS5IclmT9vn47t0OCt0zyniSXtv0vSrJT/36TvCjJSUmuS/KnJL9J8p0kj+npc69repMcRHOWF+DXPcOQ90zyrvbx88fZ10pJliY5bpqf1Tdp/kCwFvC2nu2Pe01rktcn+UWSPyS5pR0mfWiSBe36y4HnAo/qG1a9ebv+xCSXJ1k/ybeSLAVuWtY+e/a9Y5Kzk9zWfrZ7tmfOe/uc2NbQ/9p7bTvJzjRneQH+s6fOE5fz/ldMskeS89o6rm+HiT9pov0l2S7JaW3/q5P8W3/dkqT5xy9ySdIo+BrwbOBFNAGYJA8B/pdmeO+BwK9oAt9bgVOTLKqqK/q280lgZeArwO3AW4CDklxSVT9rt/tc4CjgHOBTwB+AvwWeDzwauGiCGr8CrAq8DHgncF3bfjbwu3bfu9AM1e71MmC19j1O11eBD9J8Tv86UackrwUOBv4f8BHgTzSf4wuBhwNLgHfQvP812vcz5vyexw8GTgJ+1u734ZOo8cXttr8I/B74e5o/FjwKeMMkXt/vZJrP9gPA/u17ArhmOa87FHgVcCywH7AmzR8Lfp7kOVX1y77+29IcW1+mOd5eArwHuKHdvyRpnjL0SpJGwdntz8f0tP0LsD6wad81rAfRBNaPATv3bWcl4BlVdUfb91vAZcDuNMENmjBzP2Drqrq257UfX1aBVfXzJGfThNjv9V/Tm+Q7wMuTrF5VS3tW7UITnL6zrO1PRlVdnuRm7v05jeflwM3Aln3XuX64Z1vfayexWnmc4dRjHgZ8oqo+NECZG9H8G5wBkGRfmve+c5KvVNUpA2yLqrosybE0offny6j1L5JsRRN4jwR26Lkm+wjgDOALwHP6XvYE4Alj/65JvkxznP0Thl5Jmtcc3ixJGgU3tT9Xhb9MTvQamrN8v0uyxtgC3AKcQjO5U78vjQVegKr6Hc2Z2w17+tzY/nzFkIeu7k8Tul8z1tAOjX4ecGhVDeu605toP6dluJFmEqwXtZ/ldPz7gP2PHQu88Jdrbvdqn75smrVM1th+PjEWeNtazgZ+ADx7bJh3j3v9IaN93QnAmkkePMP1SpJmkKFXkjQKxkLcWPhdQHOWcWuaobj9y1bAI8bZzmXjtF3fbmvMvsAvgS8BS5McneTt44SggVTViTQBe5ee5jfQzFb91elsu8+q3PM5TeSTwBXA94AlSb6d5E1JVhlwX0umMKP2+eO0ndf+XH+cdTNhPeDuCWo5t6dPr4mOHbj38SNJmmcMvZKkUfDk9ueF7c+xs5M/oQm44y0vGGc7d02w/b+c7ayq62lu/7MF8B80t0/aB7goyTOn/hYAOAB4SpKnJ7kfzfDrxb3Ds6ejPXO8Cvd8TuOqqouBx9Nc+3swzfW0BwAXJNlggF3eOoUyJ3svxIn6DePs+1TObk907Ex1e5KkEeE1vZKkUTB2dvSH7c8lNBNMrVpV/RNDTVtV3QWc2C4keTJwOvAhmqA44UuXs+mDgE/QvJ/v00we9alpFXtvb2p//nCZvYCqup3mvshHAyTZtn3du7hn9ufJBtRBPH4Zbb1nU5cCTx+n73hngwet81KaP4o8jnuuF++v5dcDblOSNE95pleSNGeSrJDk32lmbj56bIblqrqbZvbdjZNsP8FrJzOT8HivW2Oc5gtoZjhefTkv/2P7c9x+VXUdzZDif6CZPOtW4BtTqbNfmvv0vhe4imZm5GX1He89jl1n21v7H4HVhnDdb6+tkjytp5bQ1A3NZzPmImCVJBv39L0f955JurdOWP6/z5ix/by/970leSLNbNI/raolk9yWJGme80yvJGm2PK29lQ40Q3QfC7yUZujtMTRBsdcHgWcBRyY5kmbyqjva/tvSnJndeQp1HJBknXafV9Dc4ujVbU2HLOe1YzMPfybJocBtwLlVdW5Pn/1pZg7eDji4qpZ3/W2/tXs+p5VpzhZvDWwMXAK8fBLX2R6T5EaaicB+CzyU5rMq4Ot972c7YN8k/0szxPf4vlmtB3UWcHySLwJX08yW/Xzg61X1855++wPvBr6b5PM0/7bbM/7/m5xHMxv1W5PcSjMK4NqqOn68Aqrq2PaY2YEm1P+Ae25ZdBvw9mm8P0nSPGPolSTNlh3b5W6aM3dX0twD9rCq+lF/56q6McmzaILRq2jC053t637K1CeH+jpNANyJZsKsm2hC1fZV9e1lvbCqfpZkD2A3mmtkV6S5dVJv6D2eJpw+mqndm3cj7gmmt9AM9T6LZsj0YVX1p0lsYz+az+zNNGdHr6eZvOufquqEnn6foxlOvH37nu5Hc63zdELvUTTXHL+f5g8b19LcDupet4Sqql8neSnNpFsfb2v8Os09ci/o6/unJDvQ3Jv4czSzZJ9E81lP5DU0Z7d3Bvam+SxPAj5cVedM4/1JkuaZ9MzkL0mShiDJr4AVqurv5roWSZLu67ymV5KkIUqyJc1kSfvPdS2SJMkzvZIkDUUbdjegGdb7YGCDqrp5bquSJEle0ytJ0nB8hGYW6vOAnQy8kiSNBs/0SpIkSZI6y2t6JUmSJEmd1bnhzWussUYtXLhwrsuQJEmSJM2A008//bqqWjDZ/p0LvQsXLmTx4sVzXYYkSZIkaQYkuWKQ/g5vliRJkiR1lqFXkiRJktRZhl5JkiRJUmcZeiVJkiRJnWXolSRJkiR1lqFXkiRJktRZhl5JkiRJUmcZeiVJkiRJnWXolSRJkiR1lqFXkiRJktRZhl5JkiRJUmcZeiVJkiRJnWXolSRJkiR11opzXcB90TdO/c1clzDv/cMm6851CZIkSZLmAc/0SpIkSZI6y9ArSZIkSeosQ68kSZIkqbMMvZIkSZKkzjL0SpIkSZI6y9ArSZIkSeqsWQ29SQ5Mcm2Sc3vajkhyZrtcnuTMtn1hkj/1rPvybNYqSZIkSZr/Zvs+vQcB+wKHjDVU1avHHifZG7ixp/+lVbXRrFUnSZIkSeqUWQ29VXVykoXjrUsS4FXAlrNZkyRJkiSpu0bpmt7nANdU1cU9besl+WWSk5I8Z6IXJtk1yeIki5csWTLzlUqSJEmS5oVRCr07Aof1PL8aWLeqngq8C/hGklXHe2FV7V9Vi6pq0YIFC2ahVEmSJEnSfDASoTfJisDLgSPG2qrq9qq6vn18OnAp8Ji5qVCSJEmSNB+NROgFng9cUFVXjjUkWZBkhfbx+sCGwGVzVJ8kSZIkaR6a7VsWHQb8HHhskiuT7NKu2oF7D20G2Aw4O8lZwLeA3apq6exVK0mSJEma72Z79uYdJ2jfeZy2bwPfnumaJEmSJEndNSrDmyVJkiRJGjpDryRJkiSpswy9kiRJkqTOMvRKkiRJkjrL0CtJkiRJ6ixDryRJkiSpswy9kiRJkqTOMvRKkiRJkjrL0CtJkiRJ6ixDryRJkiSpswy9kiRJkqTOMvRKkiRJkjrL0CtJkiRJ6ixDryRJkiSpswy9kiRJkqTOMvRKkiRJkjrL0CtJkiRJ6ixDryRJkiSpswy9kiRJkqTOMvRKkiRJkjrL0CtJkiRJ6ixDryRJkiSpswy9kiRJkqTOMvRKkiRJkjrL0CtJkiRJ6ixDryRJkiSpswy9kiRJkqTOMvRKkiRJkjrL0CtJkiRJ6ixDryRJkiSpswy9kiRJkqTOMvRKkiRJkjrL0CtJkiRJ6ixDryRJkiSps2Y19CY5MMm1Sc7tadszye+SnNku2/ase3+SS5JcmOQFs1mrJEmSJGn+m+0zvQcB24zTvk9VbdQuRwMkeTywA/CE9jVfSrLCrFUqSZIkSZr3ZjX0VtXJwNJJdn8JcHhV3V5VvwYuATaeseIkSZIkSZ0zKtf07p7k7Hb482pt29rAb3v6XNm2/ZUkuyZZnGTxkiVLZrpWSZIkSdI8MQqhdz9gA2Aj4Gpg77Y94/St8TZQVftX1aKqWrRgwYKZqVKSJEmSNO/Meeitqmuq6q6quhs4gHuGMF8JPLKn6zrAVbNdnyRJkiRp/prz0JtkrZ6nLwPGZnY+CtghyUpJ1gM2BH4x2/VJkiRJkuavFWdzZ0kOAzYH1khyJfBRYPMkG9EMXb4ceDNAVf0qyZHAecCdwNuq6q7ZrFeSJEmSNL/Nauitqh3Haf7aMvp/AvjEzFUkSZIkSeqyOR/eLEmSJEnSTDH0SpIkSZI6y9ArSZIkSeosQ68kSZIkqbMMvZIkSZKkzjL0SpIkSZI6y9ArSZIkSeosQ68kSZIkqbMMvZIkSZKkzjL0SpIkSZI6y9ArSZIkSeosQ68kSZIkqbMMvZIkSZKkzjL0SpIkSZI6y9ArSZIkSeosQ68kSZIkqbMMvZIkSZKkzjL0SpIkSZI6y9ArSZIkSeosQ68kSZIkqbMMvZIkSZKkzjL0SpIkSZI6y9ArSZIkSeosQ68kSZIkqbMMvZIkSZKkzjL0SpIkSZI6y9ArSZIkSeosQ68kSZIkqbMMvZIkSZKkzjL0SpIkSZI6y9ArSZIkSeosQ68kSZIkqbMMvZIkSZKkzjL0SpIkSZI6a1ZDb5IDk1yb5Nyetn9LckGSs5N8N8lD2/aFSf6U5Mx2+fJs1ipJkiRJmv9m+0zvQcA2fW3HAk+sqicDFwHv71l3aVVt1C67zVKNkiRJkqSOmNXQW1UnA0v72o6pqjvbp6cA68xmTZIkSZKk7hq1a3rfCPxPz/P1kvwyyUlJnjPRi5LsmmRxksVLliyZ+SolSZIkSfPCyITeJB8E7gQObZuuBtatqqcC7wK+kWTV8V5bVftX1aKqWrRgwYLZKViSJEmSNPJGIvQm2QnYDnhNVRVAVd1eVde3j08HLgUeM3dVSpIkSZLmmzkPvUm2AfYA/r6qbu1pX5Bkhfbx+sCGwGVzU6UkSZIkaT5acTZ3luQwYHNgjSRXAh+lma15JeDYJACntDM1bwb8S5I7gbuA3apq6bgbliRJkiRpHLMaeqtqx3GavzZB328D357ZiiRJkiRJXTbnw5slSZIkSZophl5JkiRJUmcZeiVJkiRJnWXolSRJkiR1lqFXkiRJktRZhl5JkiRJUmcZeiVJkiRJnTVQ6E3y4JkqRJIkSZKkYRv0TO/vk3wlyVNnpBpJkiRJkoZo0ND7OWA7YHGS05K8MckDZ6AuSZIkSZKmbaDQW1UfAtYFXgksBQ4ArkryH0meOAP1SZIkSZI0ZQNPZFVVd1XVd6rqBcCGwH40IfisJD9N8tok9x92oZIkSZIkDWq6szdfD1xFc9Y3wJrAwcAlSZ45zW1LkiRJkjQtUwq9STZJciDwO+CTwM+Ap1XVo4EnAb8B9h9alZIkSZIkTcGKg3RO8hbgzTTB9iLgg8BBVXXjWJ+qOi/Jh4GfDLNQSZIkSZIGNVDoBT4P/Dfw7qo6bhn9LgY+NeWqJEmSJEkagkFD78Kqump5narqd8CHp1aSJEmSJEnDMeg1vQ9M8uzxViR5VpINhlCTJEmSJElDMWjo/Tzw8gnWvQzYZ3rlSJIkSZI0PIOG3mcAJ06w7kRgk+kUI0mSJEnSMA0aelcFbptg3R3AQ6ZXjiRJkiRJwzNo6L0M2GKCdVsAV0yvHEmSJEmShmfQ0PtfwLuSvDnJ/QGS3D/Jm4F3AocMu0BJkiRJkqZq0FsWfQbYGNgP2DfJEmCNdjvfAz493PIkSZIkSZq6gUJvVd0FvDTJ1sBWwMOA64BjquonM1CfJEmSJElTNuiZXgCq6hjgmCHXIkmSJEnSUE0p9AIkWR34m/72qrpqWhVJkiRJkjQkA4XeJKsAnwV2AB44QbcVpluUJEmSJEnDMOiZ3n2BVwEHAecAtw+7IEmSJEmShmXQ0PtC4L1V9R8zUYwkSZIkScM06H167wecPxOFSJIkSZI0bIOG3iOBF81EIZIkSZIkDdugw5t/AHwhyYOAo4Gl/R2q6uRhFCZJkiRJ0nRNJfQCrA+8CaiedWmfO3uzJEmSJGkkDBp6t5ruDpMcCGwHXFtVT2zbVgeOABYClwOvqqobkgT4PLAtcCuwc1WdMd0aJEmSJEn3DQOF3qo6bgj7PIjm1keH9LS9Dziuqj6d5H3t8z1oZovesF02AfZrf0qSJEmStFyDTmQFQJLVkmyT5DVJVmvb7j+Z17bX/PZfC/wS4OD28cHAS3vaD6nGKcBDk6w1lZolSZIkSfc9A4feJJ8CrqKZyOoQYL121Q+TfGiKdTyiqq4GaH8+vG1fG/htT78r27b+mnZNsjjJ4iVLlkyxBEmSJElS1wwUepPsAbwT+BTwLJrJq8b8N8O/nVHGaau/aqjav6oWVdWiBQsWDLkESZIkSdJ8NehEVrsCH6+qTyTpn6X5YuDRU6zjmiRrVdXV7fDla9v2K4FH9vRbh+YssyRJkiRJyzXo8OZ1gP+dYN0dwIOnWMdRwE7t452A7/e0vz6NTYEbx4ZBS5IkSZK0PIOG3quAJ0yw7kk0txtapiSHAT8HHpvkyiS7AJ8GtkpyMc1tkT7ddj8auAy4BDgAeOuA9UqSJEmS7sMGHd78LeAjSU4DFrdtlWQD4D3A15a3garacYJVzxunbwFvG7BGSZIkSZKAwc/07klz1vV/gfPbtsOBc4Ff00xwJUmSJEnSSBjoTG9V3ZJkM+B1wAtoJpq6HtiL5n66fx5+iZIkSZIkTc2gw5upqjuB/2wXSZIkSZJG1qDDmyVJkiRJmjcGOtPbzq5cy+hSVfXY6ZUkSZIkSdJwDDq8+VT+OvQ+DNgUuAk4eRhFSZIkSZI0DINOZPXa8dqTrA78CPjhMIqSJEmSJGkYhnJNb1UtpZnB+aPD2J4kSZIkScMwzImsbgXWHeL2JEmSJEmaloFvWdQvyf2AxwMfAc6fdkWSJEmSJA3JoLM3/5m/nsjqfkCAPwIvGlJdkiRJkiRN26Bnej/DX4fe24ArgB9W1Q1DqUqSJEmSpCEYdPbmD81UIZIkSZIkDdswJ7KSJEmSJGmkDHpN7/4DdK+qevOA9UiSJEmSNDSDXtP7QmAVYFXgbuAGYDWaM8Y3ATf39O2/9leSJEmSpFk16PDmV9GE29cCK1fVAmBl4HVt+yur6pHt4j17JUmSJElzatAzvfsAe1XVN8YaqurPwKFJVgc+D2wyxPokSZIkSZqyQc/0PgW4cIJ1FwJPml45kiRJkiQNz6Ch9xpg+wnWvRK4dnrlSJIkSZI0PIMOb/48sHeSNYFv0oTgR9Bc6/si4N3DLU+SJEmSpKkbKPRW1T5JbgU+DLy4Z9VVwFuqapBbGkmSJEmSNKMGPdNLVX0lyQHAo4C1gKuBK6rq7mEXJ0mSJEnSdAwcegHagPvrdpEkSZIkaSQNOpEVSZ6c5Mgkv09yR5Knte3/mmTr4ZcoSZIkSdLUDBR6k/wf4FSaWxd9B1ihb1u7Da80SZIkSZKmZ9AzvZ8BjgMeB7wdSM+6xcDTh1SXJEmSJEnTNug1vU8HXlFVdydJ37rraG5fJEmSJEnSSBj0TO/twMoTrFsTuHF65UiSJEmSNDyDht6fAm9P0vu6an++EThhKFVJkiRJkjQEgw5v/ghN8P0l8E2awPvaJHsBmwIbD7c8SZIkSZKmbqAzvVX1S2Bz4A/AnjQTWb0D+Btgi6o6f8j1SZIkSZI0ZYOe6aWqTgOem+SBwBrADVV189ArkyRJkiRpmiZ9pjfJA5Jcm+TFAFV1a1X9xsArSZIkSRpVkz7TW1V3tLcpum3YRSR5LHBET9P6NNcPPxT4R2BJ2/6Bqjp62PuXJEmSJHXToLM3HwW8YthFVNWFVbVRVW1Ecy/gW4Hvtqv3GVtn4JUkSZIkDWLQa3qPAvZNcjjwPeBq7rllEQBVdfI0a3oecGlVXdGcWJYkSZIkaWoGDb1jZ19f1S69gTft8xWmWdMOwGE9z3dP8npgMfDuqrqh/wVJdgV2BVh33XWnuXtJkiRJUlcMGnq3mpEqWkkeAPw98P62aT/g4zRh+uPA3sAb+19XVfsD+wMsWrSo+tdLkiRJku6blht6k2wJ/KKq/lhVx81wPS8EzqiqawDGfrZ1HAD8YIb3L0mSJEnqkMlMZHUs8PixJ0nul+TkJBvOQD070jO0OclaPeteBpw7A/uUJEmSJHXUZIY3988mFeDZwCrDLCTJA2mGT7+5p3mvJBvRDG++vG+dJEmSJEnLNOg1vTOmqm4FHtbX9ro5KkeSJEmS1AGD3qdXkiRJkqR5Y7JnetdOsn77eIWetj/0d6yqy4ZSmSRJkiRJ0zTZ0Putcdq+N0Hf6d6nV5IkSZKkoZhM6H3DjFchSZIkSdIMWG7oraqDZ6MQSZIkSZKGzYmsJEmSJEmdZeiVJEmSJHWWoVeSJEmS1FmGXkmSJElSZxl6JUmSJEmdZeiVJEmSJHWWoVeSJEmS1FmGXkmSJElSZxl6JUmSJEmdZeiVJEmSJHWWoVeSJEmS1FmGXkmSJElSZxl6JUmSJEmdZeiVJEmSJHWWoVeSJEmS1FmGXkmSJElSZxl6JUmSJEmdZeiVJEmSJHWWoVeSJEmS1FmGXkmSJElSZxl6JUmSJEmdZeiVJEmSJHWWoVeSJEmS1FmGXkmSJElSZxl6JUmSJEmdZeiVJEmSJHWWoVeSJEmS1FmGXkmSJElSZxl6JUmSJEmdteJcF9AryeXAzcBdwJ1VtSjJ6sARwELgcuBVVXXDXNUoSZIkSZo/RvFM7xZVtVFVLWqfvw84rqo2BI5rn0uSJEmStFyjGHr7vQQ4uH18MPDSOaxFkiRJkjSPjFroLeCYJKcn2bVte0RVXQ3Q/nx4/4uS7JpkcZLFS5YsmcVyJUmSJEmjbKSu6QWeVVVXJXk4cGySCybzoqraH9gfYNGiRTWTBUqSJEmS5o+ROtNbVVe1P68FvgtsDFyTZC2A9ue1c1ehJEmSJGk+GZnQm+RBSVYZewxsDZwLHAXs1HbbCfj+3FQoSZIkSZpvRml48yOA7yaBpq5vVNWPkpwGHJlkF+A3wCvnsEZJkiRJ0jwyMqG3qi4DnjJO+/XA82a/IkmSJEnSfDcyw5slSZIkSRo2Q68kSZIkqbMMvZIkSZKkzjL0SpIkSZI6y9ArSZIkSeosQ68kSZIkqbMMvZIkSZKkzjL0SpIkSZI6y9ArSZIkSeosQ68kSZIkqbMMvZIkSZKkzjL0SpIkSZI6y9ArSZIkSeosQ68kSZIkqbMMvZIkSZKkzjL0SpIkSZI6y9ArSZIkSeosQ68kSZIkqbMMvZIkSZKkzjL0SpIkSZI6y9ArSZIkSeosQ68kSZIkqbMMvZIkSZKkzjL0SpIkSZI6y9ArSZIkSeosQ68kSZIkqbMMvZIkSZKkzjL0SpIkSZI6y9ArSZIkSeosQ68kSZIkqbMMvZIkSZKkzjL0SpIkSZI6y9ArSZIkSeqskQi9SR6Z5IQk5yf5VZJ/btv3TPK7JGe2y7ZzXaskSZIkaf5Yca4LaN0JvLuqzkiyCnB6kmPbdftU1b/PYW2SJEmSpHlqJEJvVV0NXN0+vjnJ+cDac1uVJEmSJGm+G4nhzb2SLASeCpzaNu2e5OwkByZZbc4KkyRJkiTNOyMVepM8GPg28I6qugnYD9gA2IjmTPDeE7xu1ySLkyxesmTJrNUrSZIkSRptIxN6k9yfJvAeWlXfAaiqa6rqrqq6GzgA2Hi811bV/lW1qKoWLViwYPaKliRJkiSNtJEIvUkCfA04v6o+29O+Vk+3lwHnznZtkiRJkqT5ayQmsgKeBbwOOCfJmW3bB4Adk2wEFHA58Oa5KU+SJEmSNB+NROitqp8CGWfV0bNdiyRJkiSpO0ZieLMkSZIkSTPB0CtJkiRJ6ixDryRJkiSpswy9kiRJkqTOMvRKkiRJkjrL0CtJkiRJ6ixDryRJkiSpswy9kiRJkqTOMvRKkiRJkjrL0CtJkiRJ6ixDryRJkiSpswy9kiRJkqTOMvRKkiRJkjrL0CtJkiRJ6ixDryRJkiSpswy9kiRJkqTOMvRKkiRJkjrL0CtJkiRJ6ixDryRJkiSpswy9kiRJkqTOMvRKkiRJkjrL0CtJkiRJ6ixDryRJkiSpswy9kiRJkqTOMvRKkiRJkjrL0CtJkiRJ6ixDryRJkiSpswy9kiRJkqTOMvRKkiRJkjrL0CtJkiRJ6ixDryRJkiSpswy9kiRJkqTOMvRKkiRJkjrL0CtJkiRJ6qwV57qAyUiyDfB5YAXgq1X16TkuSZL4xqm/mesSOuEfNll3rkuQJEkdNvJnepOsAHwReCHweGDHJI+f26okSZIkSfPBfDjTuzFwSVVdBpDkcOAlwHlzWpXmlGfYps+za5J0D1yiTsQAABBOSURBVH+vTJ+/VzQK/G95OLr23/PIn+kF1gZ+2/P8yrZNkiRJkqRlmg9nejNOW92rQ7IrsGv79I9JLpzxqqZuDeC6uS5Ces34zR6fmnUTHIvj8fjUKLvPH58D/Les2XefPz41mFn+73kqx+ejBuk8H0LvlcAje56vA1zV26Gq9gf2n82ipirJ4qpaNNd1SOPx+NQo8/jUKPP41Cjz+NQom43jcz4Mbz4N2DDJekkeAOwAHDXHNUmSJEmS5oGRP9NbVXcm2R34Mc0tiw6sql/NcVmSJEmSpHlg5EMvQFUdDRw913UMybwYhq37LI9PjTKPT40yj0+NMo9PjbIZPz5TVcvvJUmSJEnSPDQfrumVJEmSJGlKDL2zKMk2SS5MckmS9811PequJJcnOSfJmUkWt22rJzk2ycXtz9Xa9iT5Qntcnp3kaT3b2antf3GSnXran95u/5L2tePdWkwCIMmBSa5Ncm5P24wfjxPtQ+o1wfG5Z5Lftd+hZybZtmfd+9tj7cIkL+hpH/d3fDsR56ntcXhEOyknSVZqn1/Srl84O+9Y80mSRyY5Icn5SX6V5J/bdr9DNeeWcXyO3ndoVbnMwkIzCdelwPrAA4CzgMfPdV0u3VyAy4E1+tr2At7XPn4f8Jn28bbA/9DcE3tT4NS2fXXgsvbnau3j1dp1vwCe2b7mf4AXzvV7dhndBdgMeBpwbk/bjB+PE+3DxaV3meD43BN4zzh9H9/+/l4JWK/9vb7Csn7HA0cCO7SPvwy8pX38VuDL7eMdgCPm+rNwGb0FWAt4Wvt4FeCi9jj0O9RlzpdlHJ8j9x3qmd7ZszFwSVVdVlV3AIcDL5njmnTf8hLg4PbxwcBLe9oPqcYpwEOTrAW8ADi2qpZW1Q3AscA27bpVq+rn1XzTHNKzLemvVNXJwNK+5tk4Hifah/QXExyfE3kJcHhV3V5VvwYuofn9Pu7v+PaM2ZbAt9rX9x/rY8fnt4DnOWpG/arq6qo6o318M3A+sDZ+h2oELOP4nMicfYcaemfP2sBve55fybIPCmk6CjgmyelJdm3bHlFVV0PzJQU8vG2f6NhcVvuV47RLg5iN43GifUiTsXs7PPTAnmGdgx6fDwP+UFV39rXfa1vt+hvb/tK42uGbTwVOxe9QjZi+4xNG7DvU0Dt7xvvLg1Nna6Y8q6qeBrwQeFuSzZbRd6Jjc9B2aRg8HjUK9gM2ADYCrgb2btuHeXx67GrSkjwY+Dbwjqq6aVldx2nzO1Qzapzjc+S+Qw29s+dK4JE9z9cBrpqjWtRxVXVV+/Na4Ls0w0auaYcx0f68tu0+0bG5rPZ1xmmXBjEbx+NE+5CWqaquqaq7qupu4ACa71AY/Pi8jmZ46Yp97ffaVrv+IUx+mLXuQ5LcnyZQHFpV32mb/Q7VSBjv+BzF71BD7+w5DdiwnYHsATQXXB81xzWpg5I8KMkqY4+BrYFzaY63sdkadwK+3z4+Cnh9O+PjpsCN7TCmHwNbJ1mtHZayNfDjdt3NSTZtr514fc+2pMmajeNxon1IyzT2P/qtl9F8h0JzTO3Qzhq6HrAhzSRA4/6Ob6+RPAHYvn19/7E+dnxuDxzf9pf+ov1e+xpwflV9tmeV36GacxMdnyP5HToTM3m5TDjD2bY0s5pdCnxwrutx6eZCM/PdWe3yq7FjjeY6h+OAi9ufq7ftAb7YHpfnAIt6tvVGmkkGLgHe0NO+qP0CuxTYF8hcv2+X0V2Aw2iGN/2Z5i+zu8zG8TjRPlxcepcJjs+vt8ff2e3/WK3V0/+D7bF2IT0z10/0O779Tv5Fe9x+E1ipbf+b9vkl7fr15/qzcBm9BXg2zZDNs4Ez22Vbv0NdRmFZxvE5ct+hYwe1JEmSJEmd4/BmSZIkSVJnGXolSZIkSZ1l6JUkSZIkdZahV5IkSZLUWYZeSZIkSVJnGXolSX+R5JlJjkxyVZI7klyf5NgkOyVZYY5rW5hkzyTrD3m7ayY5KsnSJJXkHeP0ObFdt7xl4TBrm6ok2yd5+1zXMUztfR3fnuTnSW5McnuSy5J8NcmTe/qdkuRHM1TDe5L8/UxsW5I0c1ac6wIkSaOhDXufBY4H9gCuAFYDtgb2A/7APTeFnwsLgY8CPwUuG+J2PwI8F9iZ5n6tl4/T563Aqj3PPww8A+gPQFcPsa7p2J7m3ptfmOtChiHJqsCPgSfTHIsfB24BHgO8rl231iyU8h7gBzT3nZQkzROGXkkSSTajCbz7VlX/GcLvJ/ks8KDZr2xWPA44q6q+O1GHqjqv93mSJcAdVXXKTBd3X5Fkpaq6fYLVXwI2Ap5dVaf3tJ8EHJDkZTNe4AxZzvuWJA2Bw5slSQDvA5YC7x1vZVVdWlVnjz1PsnGSnyT5Y5JbkhyXZOPe17RDgk/s31aSy5Mc1PN853Zo8KZJDk1yUzu8+gtJ/qbtszlwQvuSY3uGE28+0RtK451JLmyHal+dZN/2rOHYcOkCNgeeM4whyknWabexfU/bK9u2r/a0PSTJXUne2NP2iCQHtHXenuS8JG8YZx+PTnJ4kuuS3Jbk9CQv6ll/OPBqYIOe93RBz36/lOS37T6uSXJMkkcv5339vh1G/NZ2SPFtSU5L8pxx+j6//bf/Y7v8MMnj+vqc0h4/L09yVpLbgTf2b6vtuxDYEfhiX+D9i2X9wSLJbu1nsGZf+6eT3NbX9p4kFyT5U5rh7r9Ist3YZwA8Atil53P9cs9rn57kB0n+0L7+5CTP7Nv+4UkuSbJZ+xn8CfiXiWqXJA2HoVeS7uPSXKu7OXBMVd22nO6010+eRDP0eWfg9TRDf09K8pRplPJ14FLg5TRDWN8GvL9dd0b7HODtwDPb5YxlbO8TNGevjwVeDOzV1vvDJPejGYr8TOBs4Jc925zyEOWquhK4GNiyp3lL4E99bZvT/A4+ASDJasDPgecBHwK2A44BvpbkH8delOZ65lOBv6P5HF4CnAccleQFbbcPAT8Brux5T69u1+3bvuYjwFbAW9rX9w7dnsgL2v57AP/Qtv04yXo99b2cZqjxdW2f1wELgJOT9A8/fiLwbzT/RtsA/2+C/W5J81nN6JDiJLsAnwQOAbalqf17wOptl21p/jB0FPd8rp9pX7spzbD7BwG70AwvvwU4PsmT+na1Bs2xfgjwQuBbM/amJEmAw5slSc3/hK9Mcw3vZHwEuB14XlX9ASDJsTTXwn6UJrROxTeq6qPt458k2YTmDN9Hq+qmJGNDjM9f3rDiJKsD7wIOrqrd2+YfpxmW/HVgu6o6Cjglyc3AnUMcqnwCsEXP8y1oQvy7kjyqqq5o266oql+3fd5DcxbxCVV1edt2bJKHAR9L8rWqupvmWtbbgedW1Y097+tRwMeAH1fVJUmuB24f5z09k+Yz+c+etu9M8n0tAJ5RVb8HSHICzTHzAeAf2z8kfL6tofdM90k012D/M82Igt7tbVFV5y9nv49sf072+JyqZwKLq+qTPW0/HHtQVWck+TOwZJzPdW/gImCrqroTIMkxwAXAB4Edevo+BHh1Vf14Bt6DJGkcnumVJA1qM+AHY4EXoKpuojkD9txpbPeHfc/PAdad4rY2BVYC/quv/XDgTqZX5/KcAPxdkrWS/C3wWOAA4Dfcc7Z3S5oJw8ZsQ3Om8MokK44t3DNB06N7+v03cEtfv2OAZ6QdDr4MpwG7JtkjydPaoDpZJ48FXoCquqGtb2wI7xOAdYD/6qvtpna/m/Vt78JJBN7ZdBqwSZJ9kmyZZOXJvCjNcPlnAke0z8fed9H8G/e/71sNvJI0uwy9kqTraYbfPmqS/Vdn/CHAv6cZ8jxVS/ue304TXKdibEjqvepsz8Jd37N+Joxde7xFu1xdVRe07VskWYNmaO8JPa95OM0s2X/uW77ern9YOwx9dWDXcfp9nOZ3+vI+/zcDB7Y/TweuSfJvkwjLANdM0LZ2z3sAOHSc+p4PPKzvtZMdRv7b9udkj8+pOoBmyPhzaIaHL03yzSTrLOd1C4DQDKfvf99v4q/f9++RJM0qhzdL0n1cVd2ZZsKprTK5mWSXAmuO074m9w6utzH+taIzGTjHjNWxJvCrscb2DNzDaILvjKiqa9qh2FvShKGxcHs8zTWjW/S109ZzCfB/J9jsBVV1V5Ibac6I7zNBv+uWU9tNNJOVvbe9FvdVNGHtVpqh6cvyiAnaftc+HvtM3w2cPE7f/uvFazn7G3M8cDfNddnjbXd5xvb7gL72e4XRdvj4F4EvtsPKt6EZtnwoyx4ZMHas7U0zkqBf//uc7PuWJA2JoVeSBPBp4ESaiYX6b1lEG5BWaWdwPgl4UZJVqurmdv0qNKHkxJ6XXQG8IskDquqOtt9mwCpTrHEsjE9m2Okpbf8dgON62l9N87vvpCnWMFkn0ExSNHYGEJrwtjbNWdaL20mvxvyIZpKty6qq/4w3ff2eApyznD9O3M5yPqf2euLPJNmJ5szz8jwnyZo91/SuRjO51RHt+nOAq4DHVdVnJ7G9Samqy5McBuye5PDxZnBO8tKq+t4Emxi7FviJNEPMSbISzaRhE+3zeuDQJM/inknAYJzPtapuSHIqzT2E/29VGWolacQYeiVJVNXJSd4FfLa9vcxBNAFhNZpw8Caa2XjPphlKux1wXJLP0Jy52gN4IPe+/crhNENxD0xzi6L1aCaXupGpuYjmetw3JllKE0AuHAvefe9naZp7C78/yS3A0TT34/1Xmmtn+68fHrbjuWe26RPamq5McgnN57l/X/+9aGb8/WmSz9G811Xamjepqle0/T5AM3vziUm+xD3/Rk8G/raqdmv7nQe8vp2R+Gya60h/lWQxcCTN2e9baIYd/x3NBFTLcx3N5Fr/AtxFM7P2irShvj0TvTvwzSQPBL5Nc/Z3TeBZwEVVte8k9jOetwEb0MwQvh/N53tr2/a69j1MFHp/RjNEep827N4N/BN9l3i1x+gSmj+YLGm3uQPNdctjzqMZor4tcC1wbVX9BnhHW9PR7XZ+TzPseRHw56r68BTftyRpCAy9kiQAqupzSX4BvBP4d5pZnW8GFtOcnfzvtt/Zae6P+wngYJqzmafQzCh8Vs/2TkiyG83MxK+guS3Qa2nC0FTqu74NVXvQnKldgWao8IkTvOSDNOFlN+CtNAHsEOD97VDWmXQizR8DflNVl/W0H08zKVXv0OaxkL4JzRDjD9JMXnUDzey/R/T0uyzJImBPmtvlrEETRs8Bemdk3o8mcO1NM1vwhTQh7mSaP16sR/P5XQrsXlVfmcR7+jHNLaL2Av623ecLemabpqq+m2QLmnD+NZqzolfT3I6pf1KxSauqG9tjbre2/t1orvf+Lc2Z/N2X8do7kryY5nZNX6f5vPamOev+jp6uP6W5/dbONH9wuKp9Dx/r6fNe4Ms0txlaGfgKsFtVndL++32k3c+qNKF4Mc2/hSRpDsVROJIkaVmS/J5mxu43zXUtkiQNytmbJUmSJEmdZeiVJEmSJHWWw5slSZIkSZ3lmV5JkiRJUmcZeiVJkiRJnWXolSRJkiR1lqFXkiRJktRZhl5JkiRJUmcZeiVJkiRJnfX/AVI5k0lGPZHcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import seaborn as sns\n",
    "\n",
    "for col in ['tweetcount']: #,'usercount']:\n",
    "    plt.figure(figsize=(16, 6))\n",
    "    sns.distplot(df2[col], kde=False, bins=10)\n",
    "    plt.title('Density Distribution', fontsize=18)\n",
    "    plt.xlabel('Count of Tweets per Cluster', fontsize=16)\n",
    "    plt.ylabel('Frequency', fontsize=16)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loop through labels to determin the convex hull (approximate perimeter) around each cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for kmeans clusters\n",
    "df0['geom'] = \"\"\n",
    "for cluster in df0.labels:\n",
    "    df3 = df[df.mbk_means_labels == cluster]\n",
    "    points = [[float(lng), float(lat)] for lng, lat in zip(df3.lng, df3.lat)]\n",
    "    multi_point = MultiPoint(points)\n",
    "    convex = multi_point.convex_hull\n",
    "    hull = str(convex)\n",
    "    df0.loc[df0['labels'] == cluster, ['geom']] = hull"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for dbscan clusters\n",
    "df2['geom'] = \"\"\n",
    "for cluster in df2.label:\n",
    "    df3 = df[df.labels == cluster]\n",
    "    points = [[float(lng), float(lat)] for lng, lat in zip(df3.lng, df3.lat)]\n",
    "    multi_point = MultiPoint(points)\n",
    "    convex = multi_point.convex_hull\n",
    "    hull = str(convex)\n",
    "    df2.loc[df2['label'] == cluster, ['geom']] = hull"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Write the cluster dataframe to a psql/postgis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "df0.to_sql('mbkmclusters', engine, if_exists='replace')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "df2.to_sql('clusters', engine, if_exists='replace')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Alter the psql/postgis clusters table to have spatial geometry "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "cur.execute(\"SELECT AddGeometryColumn ('mbkmclusters', 'geometry', 4326, 'POLYGON', 2);\")\n",
    "cur.execute(\"UPDATE mbkmclusters SET geometry = ST_SetSRID(ST_GeomFromText(geom), 4326);\")\n",
    "cur.execute(\"ALTER TABLE mbkmclusters DROP COLUMN geom;\")\n",
    "conn.commit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "cur.execute(\"SELECT AddGeometryColumn ('clusters', 'geometry', 4326, 'POLYGON', 2);\")\n",
    "cur.execute(\"UPDATE clusters SET geometry = ST_SetSRID(ST_GeomFromText(geom), 4326);\")\n",
    "cur.execute(\"ALTER TABLE clusters DROP COLUMN geom;\")\n",
    "conn.commit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "conn.rollback() # use this command if error is given during query"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load shapefile of counties"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1) copy shapefile of counties\n",
    "\n",
    "2) change directory to location of shapepile\n",
    "\n",
    "3) in terminal: sudo -su postgres\n",
    "\n",
    "4) in terminal: shp2pgsql -I -W \"latin1\" -s 4326 tl_2010_06_county10/tl_2010_06_county10.shp tl_2010_06_county10 | psql -d tweets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['gid', 'statefp10', 'countyfp10', 'countyns10', 'geoid10', 'name10', 'namelsad10', 'lsad10', 'classfp10', 'mtfcc10', 'csafp10', 'cbsafp10', 'metdivfp10', 'funcstat10', 'aland10', 'awater10', 'intptlat10', 'intptlon10', 'geom']\n"
     ]
    }
   ],
   "source": [
    "cur.execute(\"Select * FROM tl_2010_06_county10;\")\n",
    "colnames = [desc[0] for desc in cur.description]\n",
    "print (colnames)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['index', 'geoid', 'sumlev', 'state', 'county', 'cbsa', 'csa', 'necta', 'cnecta', 'name', 'pop100', 'hu100', 'pop1002000', 'hu1002000', 'p001001', 'p0010012000']\n"
     ]
    }
   ],
   "source": [
    "cur.execute(\"Select * FROM ca_census_data;\")\n",
    "colnames = [desc[0] for desc in cur.description]\n",
    "print (colnames)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['index', 'label', 'tweetcount', 'usercount', 'geometry']\n"
     ]
    }
   ],
   "source": [
    "cur.execute(\"Select * FROM clusters;\")\n",
    "colnames = [desc[0] for desc in cur.description]\n",
    "print (colnames)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Copy census population data per county into posgis table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>geoid</th>\n",
       "      <th>sumlev</th>\n",
       "      <th>state</th>\n",
       "      <th>county</th>\n",
       "      <th>cbsa</th>\n",
       "      <th>csa</th>\n",
       "      <th>necta</th>\n",
       "      <th>cnecta</th>\n",
       "      <th>name</th>\n",
       "      <th>pop100</th>\n",
       "      <th>hu100</th>\n",
       "      <th>pop1002000</th>\n",
       "      <th>hu1002000</th>\n",
       "      <th>p001001</th>\n",
       "      <th>p0010012000</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>6001</td>\n",
       "      <td>50</td>\n",
       "      <td>6</td>\n",
       "      <td>1</td>\n",
       "      <td>41860</td>\n",
       "      <td>488</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Alameda County</td>\n",
       "      <td>1510271</td>\n",
       "      <td>582549</td>\n",
       "      <td>1443741</td>\n",
       "      <td>540183</td>\n",
       "      <td>1510271</td>\n",
       "      <td>1443741</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>6003</td>\n",
       "      <td>50</td>\n",
       "      <td>6</td>\n",
       "      <td>3</td>\n",
       "      <td>99999</td>\n",
       "      <td>999</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Alpine County</td>\n",
       "      <td>1175</td>\n",
       "      <td>1760</td>\n",
       "      <td>1208</td>\n",
       "      <td>1514</td>\n",
       "      <td>1175</td>\n",
       "      <td>1208</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>6005</td>\n",
       "      <td>50</td>\n",
       "      <td>6</td>\n",
       "      <td>5</td>\n",
       "      <td>99999</td>\n",
       "      <td>999</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Amador County</td>\n",
       "      <td>38091</td>\n",
       "      <td>18032</td>\n",
       "      <td>35100</td>\n",
       "      <td>15035</td>\n",
       "      <td>38091</td>\n",
       "      <td>35100</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>6007</td>\n",
       "      <td>50</td>\n",
       "      <td>6</td>\n",
       "      <td>7</td>\n",
       "      <td>17020</td>\n",
       "      <td>999</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Butte County</td>\n",
       "      <td>220000</td>\n",
       "      <td>95835</td>\n",
       "      <td>203171</td>\n",
       "      <td>85523</td>\n",
       "      <td>220000</td>\n",
       "      <td>203171</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>6009</td>\n",
       "      <td>50</td>\n",
       "      <td>6</td>\n",
       "      <td>9</td>\n",
       "      <td>99999</td>\n",
       "      <td>999</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Calaveras County</td>\n",
       "      <td>45578</td>\n",
       "      <td>27925</td>\n",
       "      <td>40554</td>\n",
       "      <td>22946</td>\n",
       "      <td>45578</td>\n",
       "      <td>40554</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   geoid  sumlev  state  county   cbsa  csa  necta  cnecta              name  \\\n",
       "0   6001      50      6       1  41860  488    NaN     NaN    Alameda County   \n",
       "1   6003      50      6       3  99999  999    NaN     NaN     Alpine County   \n",
       "2   6005      50      6       5  99999  999    NaN     NaN     Amador County   \n",
       "3   6007      50      6       7  17020  999    NaN     NaN      Butte County   \n",
       "4   6009      50      6       9  99999  999    NaN     NaN  Calaveras County   \n",
       "\n",
       "    pop100   hu100  pop1002000  hu1002000  p001001  p0010012000  \n",
       "0  1510271  582549     1443741     540183  1510271      1443741  \n",
       "1     1175    1760        1208       1514     1175         1208  \n",
       "2    38091   18032       35100      15035    38091        35100  \n",
       "3   220000   95835      203171      85523   220000       203171  \n",
       "4    45578   27925       40554      22946    45578        40554  "
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ca_census_data = pd.read_csv('all_050_in_06.P1.csv')\n",
    "colnames = [x.lower() for x in ca_census_data.columns.values] # convert column names to lowercase\n",
    "ca_census_data.columns = colnames\n",
    "ca_census_data.to_sql('ca_census_data', engine, if_exists='replace')\n",
    "ca_census_data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Population per county"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>county</th>\n",
       "      <th>population</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Alpine County</td>\n",
       "      <td>1175</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Sierra County</td>\n",
       "      <td>3240</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Modoc County</td>\n",
       "      <td>9686</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Trinity County</td>\n",
       "      <td>13786</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Mono County</td>\n",
       "      <td>14202</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           county  population\n",
       "0   Alpine County        1175\n",
       "1   Sierra County        3240\n",
       "2    Modoc County        9686\n",
       "3  Trinity County       13786\n",
       "4     Mono County       14202"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sql = '''\n",
    "SELECT name as county, pop100 AS population\n",
    "FROM ca_census_data\n",
    "ORDER BY pop100;'''\n",
    "\n",
    "x = pd.read_sql_query(sql, con=conn)\n",
    "x.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Number of tweets per county"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>county</th>\n",
       "      <th>tweetcount</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Alameda County</td>\n",
       "      <td>169909</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Alpine County</td>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Amador County</td>\n",
       "      <td>562</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Butte County</td>\n",
       "      <td>3659</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Calaveras County</td>\n",
       "      <td>72</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             county  tweetcount\n",
       "0    Alameda County      169909\n",
       "1     Alpine County           7\n",
       "2     Amador County         562\n",
       "3      Butte County        3659\n",
       "4  Calaveras County          72"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sql = \"\"\"\n",
    "SELECT t1.namelsad10 AS county, COUNT(T2.id) AS tweetcount\n",
    "FROM tl_2010_06_county10 T1,tweetpoints T2\n",
    "WHERE ST_Contains(T1.geom,T2.location)\n",
    "GROUP BY T1.namelsad10 \n",
    "ORDER BY T1.namelsad10;\n",
    "\"\"\"\n",
    "\n",
    "y = pd.read_sql_query(sql, con=conn)\n",
    "y.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tweets per capita in each county"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create postgis table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "sql = \"\"\"\n",
    "DROP TABLE IF EXISTS tweetspercapita;\n",
    "CREATE TABLE tweetspercapita AS\n",
    "(SELECT t3.namelsad10 AS county, t3.x AS tweetcount, t4.pop100 AS population, ROUND(t3.x*100.0/t4.pop100,2) AS tweetspercapita, t5.geom\n",
    "FROM\n",
    "    (SELECT t1.namelsad10, CAST(COUNT(t2.id) AS NUMERIC) AS x\n",
    "    FROM tl_2010_06_county10 t1,tweetpoints t2\n",
    "    WHERE ST_Contains(t1.geom,t2.location)\n",
    "    GROUP BY t1.namelsad10 \n",
    "    ORDER BY t1.namelsad10) AS t3\n",
    "    LEFT JOIN\n",
    "        (SELECT name, pop100\n",
    "        FROM ca_census_data\n",
    "        ORDER BY pop100) AS t4\n",
    "    ON t3.namelsad10 = t4.name\n",
    "    LEFT JOIN\n",
    "        (SELECT geom, namelsad10\n",
    "        FROM tl_2010_06_county10) AS t5\n",
    "    ON t4.name = t5.namelsad10);\"\"\"\n",
    "cur.execute(sql)\n",
    "conn.commit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create dataframe from query"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "sql = \"\"\"\n",
    "SELECT t3.namelsad10 AS county, t3.x AS tweetcount, t4.pop100 AS population, ROUND(t3.x*100.0/t4.pop100,2) AS tweetspercapita, t5.geom\n",
    "FROM\n",
    "    (SELECT t1.namelsad10, CAST(COUNT(t2.id) AS NUMERIC) AS x\n",
    "    FROM tl_2010_06_county10 t1,tweetpoints t2\n",
    "    WHERE ST_Contains(t1.geom,t2.location)\n",
    "    GROUP BY t1.namelsad10 \n",
    "    ORDER BY t1.namelsad10) AS t3\n",
    "    LEFT JOIN\n",
    "        (SELECT name, pop100\n",
    "        FROM ca_census_data\n",
    "        ORDER BY pop100) AS t4\n",
    "    ON t3.namelsad10 = t4.name\n",
    "    LEFT JOIN\n",
    "        (SELECT geom, namelsad10\n",
    "        FROM tl_2010_06_county10) AS t5\n",
    "    ON t4.name = t5.namelsad10;\"\"\"\n",
    "\n",
    "df4 = pd.read_sql_query(sql, con=conn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Write the dataframe to a geojson file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_counties = gpd.read_file(\"tl_2010_06_county10/tl_2010_06_county10.shp\")\n",
    "geometry = gdf_counties.geometry\n",
    "gdf = gpd.GeoDataFrame(df4, geometry=geometry)\n",
    "gdf.to_file('tweetspercapita.geojson', driver='GeoJSON')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chloropleth map of the number of tweets per capita by county"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcwAAAEvCAYAAADIC/GqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd3hcxbnA4d+cbeplJVnFkiy5SLIt9wY4gDHNtEAIKUACJHAJpEBCbgikQG4KgRSS3NCSAAlJCIQQiLlAAAN2DDa2sXFTsSVZvfdets39Y1e2ZO9qd7UraWXN+zz7WHt2zpxZW95vZ87MN0JKiaIoiqIoY9OmugGKoiiKMh2ogKkoiqIoPlABU1EURVF8oAKmoiiKovhABUxFURRF8YEKmIqiKIriA/1kXiwxMVFmZWVN5iUVRVG82rdvX6uUMmmq2xGqLj4vUra12yf9uvsODb0ppdw06Rf2YFIDZlZWFnv37p3MSyqKonglhKia6jaEsrZ2O3vezJz06+pSSxMn/aJjmNSAqSiKokw/EnDgmOpmTDkVMBVFURQvJHapAqYKmIqiKMqYnD1MlUZVBUxFURTFKzUkqwKmoiiK4oVEYlcbdaiAqSiKoninhmRVwFQURVG8kIBdBUzfA6YQQgfsBeqklJcLIZ4FVgNWYA/wJSmldWKaqSiKokwl1cP0LzXenUDxiOfPAnnAEiAcuCWI7VIURVFChATsUk76I9T4FDCFEOnAZcCTw8eklK9LF5w9zPSJaaKiKIoy1RxT8Ag1vg7J/hq4G4g++QUhhAH4PM4e6IRqrGym6IPS4891Bj2aTkOnEwhNQ4hTGodOr8NutWGz2tHpdegNOjqau7AO2ZBGQfeSOKwa2KWDcL0BvaZhsdsZsFnpt1oZsFrptgxhcziINpqINhrRac7vGQKBEGBzOJBS4nB9I3IADocDh5RYHHYsdudDSun8puaQSCQGTUOnaWhCIBCuYzqMOh1612s6IRDOv2c0IdAJgU7TMGg6NCFc5QQGTeeqCzShEVPdj+iyIARIh8RhdyA059+H0ASapmEMN2IwOv8OAaSUICVSOtsPoA3/vQqBpjn/gh0OZ7mRf89CiOPlHHYHdqsdh8PhLCals25gxflL0DSV819RphOJVPcw8SFgCiEuB5qllPuEEBvcFHkM2C6lfM/D+bcCtwJkZgaWizApPYEXHv41FQU1AdUzLH5dBrs/mRCUukLN2dt6aXi5YKqbcYo3LM+rPXIUZbqRYFfx0qePrvXAx4UQlcDzwEYhxF8BhBD3A0nAXZ5OllL+Xkq5Wkq5OikpsM0AdHodOn3wPm2NVy8IWl2hpum1I1PdBEVRThPOTD9qSNZr9JFS3iulTJdSZgGfBd6VUn5OCHELcDFwrZQTl2Swv2eAzY9v4c0/b2ewf4jVFy0LqD59nIn0/z4L+6Pn8oGuM0itDEEnD08riqKMm8A+BQ+vrRIiQwixVQhRLIQoFELc6TpuFkJsEUKUuv6M93D+ja4ypUKIG71dL5B1mE8AVcAHwnnz8CUp5Q8DqM+tx775F7b8xTna+9xDmxnoHRx3XfpoI7O+eSbvmjqhL1gtDE2n3M8NETIEZ74pijI2CThC87+uDfimlPIjIUQ0sE8IsQW4CXhHSvmgEOIe4B7g2yNPFEKYgftxLo+UrnNfkVJ2eLqYXwFTSrkN2Ob6eVKSHmQvzjj+c0N5c0B1JV2Sy66IXpj8fVAVF4fdAYapboWiKP7ypcc32aSUDUCD6+ceIUQxMBu4EtjgKvYMzrj17ZNOvxjYIqVsB3AF2k3Ac56uF/LTL66+YxM3/uAaIqLDAq6r+dUjaKHa9Qqy6IXJU90Etxwh+jVVUZTpTQiRBawAdgPJrmA6HFRnuTllNjByBmmt65hHIR8whRBc9+0reXTXjzn7E2uOL4EYD3u/lVxDTBBbF7pMs0Pzfc6Q7yuKclpxpsabknuYiUKIvSMet7prnxAiCvgn8HUpZbePb8vdp9GY3+inTS7ZtLnJfO9vd1BZWMMrv3ubt57ZjtVi87sek/+nTEsiBIdPAIRag6ko05JDTslnSquUcvVYBVy5AP4JPCulfMl1uEkIkSqlbBBCpALu7ufVcmLYFpzJd7aNda1p9+mVtTiDO/73CyxcN9/vc6PmmdllbZuAVoUgXWgGTEVRpp8p7GGOSThnnD4FFEspHx7x0ivA8KzXG4HNbk5/E7hICBHvmkV7keuYR9MuYA5bsXGx3+dE3bwsRPtdwSftobiKCaQjNNulKIpnEoEdbdIfPliPM9PcRiHEAdfjUuBB4EIhRClwoes5QojVQognAVyTfX4EfOh6/HB4ApAn02ZI9mQf+8Ranvmff/p1jnyvDnlexAS1KMRoM+WrgaIok2GKhmTHJKV8H8+rzs93U34vIzYKkVI+DTzt6/WmbQ8zIyeVTTed69c5PYXNrNdOz1R4JxMhGjDVPUxFmX5CdUh2sk3bTy8hBJtu2uDXOf3VnTR/6x0Sw8InplEhRIbo8g1dALOcFUWZKgK71Cb9EWqm7ZAsQN7aeVx371UM9AwghDi+o4am05AOB0MDVoSA5po2Pnj1IwCEhCUNOra6TZSkTLRQ7fkqiuKZM5ds6AWwyTatA6YQghvv+6TXcvvePnw8YEoBfSkmGOqd6OZNKV2ESqejKErwhOIQ6WSb1gHTVzHmKFZdsAShCcLnm6lJjWMjsQBYHXZ219ZicZxe+fJEWGj+0wqVuUBRph0pRUgOkU620PxUDbIFK7N54P/u9vj6ba9t5s3ysgm59orIJOKPDlCdp6esdxJ3RwnNW5iKokxTDtXDnBkB05ur8xYFPWAaNI1lf+/EPtROfX0v0alRxN2SRufQ+HdbOR2oHqaiTD/OWbKqh6kCJrBudgbpMTHUdvuagnBsGxyJDG2pp6PiRI8y/pr57B8ac01scKm4pChK0KghWVABE4BokwlzWETAATM7MpbsYhv1z390ymv2MF1AdfstRHtyUkrVy1SUaUbNknVSARPQhOCHG87nUy8+h9XH1G06IYg1hrGsJwLRY8PRNkjH+8eob+13X766n4TccATOXz690DAIgUAgAem66WiVDkbusexLbBEjUq3bXSeH6wdISDMjNIEtWo80aAi9BsJDin7p2tzZIY+3EeEq7JDOAzrnsh2HQ57YTXZ4mchwm4crH/FcCHGibhUsFWVasodgpp/JpgKmi/1oJxsqI7BkhnM0vJ/Ggb4xy0cajCQ/WEqD1bcAW//3YvIzzAw0ttHZNPGTf8pH/Bz+4w0cGJwhSecVRQm64VyyM50KmC5VJY3UPF8EQCyw6LZlvBvT4bH8AkMMfT4Gy2FNNe0kRKv1kYqiKNORCphAY3Ub/3rqP8efC6DuiYOc8+1VdEZKTA4BCPb2Nx8faewXduJXpdKxr8Fr/WHhRuYvTKGntYuKA5UT8RYURVEmlENN+lEBEyAlM4HvPHYT//3J32CzOhMYCKD5oX0AWFzl1uQlYL82i319LRT3tqNdYOQMXSrte8YOmoMDFprKG2kqb5rAd+GZpqbMKooSALWsxEkFTBdjmAFTuBGbdcBjme4jbcS/YmTjmcnQaaEux4C0eQ+CucsyOLbrSDCbqyiKMmkkQk36YYYHTIfDgaZpSCl5+oFX6Ov2HCyHdexrOD4MKwDPdzmd5s1PpHjrocAbe5pQy0oUZXpSy0pmaMCsr2jhlWfe49DOUq77+sUcK6jjwI6SoF9H0wRVBVVBr3c6c0iJTgVMRZlWpEQlLmCGBsxZ6Wbef/UAbU1d/ORLf5yw68QmRGGRViz9Fu+FJ1BsoxXiprQJx6kUt4oyHQmVS5ZpvIF0IPQGHfc8diM6/cS9/fBIE7EROnraeibsGr4S1tAJU1KFTEWZdiSoDaSZoT1MgMVr5pK5IIWK4voJqT9z/qzQuXcZQjHKIUOoMYqi+EzNkp3BAVMIQVpW4oQEzLBwI1oI7a8ZSgMpKl4qyvQjETjULNmZ/ZUhf928Cak3M9tM4X8KJ6Tu8QilGOWQ/mVHUhQlNNjRJv0RamZsDxMgPDIs6HXOyUmh6Vhj0OsNSAj93tlVF1NRph2JyvQDMzxgTgRTmJ6OBm+rMyeZFjpDKeoepqJMRwJ7SN3ccRJCPA1cDjRLKfNdx/4O5LqKxAGdUsrlbs6tBHoAO2CTUq72dr0ZHTAnYjlgdVkzOWfkULIr+Os6xyuUvhjafNw+TVGU0BHCPcw/AY8Afx4+IKX8zPDPQohfAl1jnH+elLLV14vN6IA5EZ0doQlM0RHBrzgQIZQowK7uYSrKtBSKPUwp5XYhRJa714QzpdingY3But6MDpjBjiN5yzNpKavn8JYDwa04UCH0e66GZBVl+pFSTFUPM1EIsXfE899LKX/v47lnA01SylIPr0vgLSGEBH7nS70zOmAGS3xSNGnp8ZR/VEZfx9gbT08FqQ+hiKkoyrQ0RYkEWn25t+jBtcBzY7y+XkpZL4SYBWwRQhyRUm4fq8KQHJSeTrJyU3B093D4nYOTHyyFQPgwoccaoZuExvhGdTAVRZloQgg9cDXwd09lpJT1rj+bgZeBtd7qndE9zPF+eOsNOtLnJmHSQdneMqyD1qC2S28ykJA9i5jUeKTNQV97D/FzZjEYFoZDSqSUhBl1HKtsJyLCSEpSFDohqG/tRdMEKWEGGsucS1tsNgc8XcVZUiIlSId0DtFK5/Co4dOZ7IjvDGr7FUU5vUiYbrlkLwCOSClr3b0ohIgENCllj+vni4Afeqt0RgfM8VqQP5vCdw76XD4yMZrMNTn0NnXQVtXCUE8/NosdIQTx6Qkk5aRhH7ISnmKm3wHlla00OCQNHa6k7YZw6up7cM6AHs1isdHZ2T/qWOq8JLraT/R2x9q2bFboJCRSFCVkiZDM7SqEeA7YgPNeZy1wv5TyKeCznDQcK4RIA56UUl4KJAMvu7Ya1AN/k1K+4e16KmD6afGqOZTsLD7+3BhpImPVPMLio5EmI3YklpoWKnaVkDQvhYSV86ms7aCwoQfQQ0YqAHq9hsMhaXNI2toGnZWVtQSljVadRu76HLoaO2mqaHH2Kj2Z5HRXuhBaE6ooim+cy0pC7/+ulPJaD8dvcnOsHrjU9XM5sMzf683wgOn/mGz1wQpSFmcSOSsOu8XKUGQEpRWt0D06WUHM2Utp6R2kqaTJbT0228Qtrzg64pr5GxdT/HbBuAZTMqzhLK020R0HDg32xHZh1QK7CSmm17COoiguoZiqbrLN8IDp/4f33ItXUVTWzEDt8FrYfrfluscYBp1MAi/v0sP/gTO64tC9VEd5Y/fxY2tmRRM/38zeswRNhqFxtUdt76Uo049Kvu40o78yaDr/3/7B1/ajq29lwYLkCWhR8AkvmXWs/6rhnLZ4TA4N4Rq6vaQyjp7fHaV9RLAE6GjuoXxnFYkPV3H5dhNxdsOEtVtRlNDiQJv0R6iZ4T3M8RnsGSJ8AjefDob0DDPm6DBq9leOWa6juQd+X8RcYO5ZWXQtCafv/UYcds+B1uGQHNtRyZqOFMo3RXHM5PtyGp0I7b83RVFOJSXYVQ9zhgfMABYFdvVbgtiQ4EtKjKDyg2P0tPsezMp3VsJO369RW9TI7J54Yq6axf6obu8nAFoIpelTFMV3akjWj4AphNABe4E6KeXlQohs4HnADHwEfF5KGdpR5CQygIAZYwidZAAjGYw6slamsS97kNzoOfS9UIxjrFmyAWqp6UD/RDdXnDeX1nk67BocieyjW2dzW179l1OU6cd5D1ONDvnTw7wTKAZiXM8fAn4lpXxeCPEEcDPweJDbF7Js/eOb9DJRFlyYxf6kLvqFjSKcy1OOZcAFucnUFU/s/pw2q52yt06ka1y3KIWu8+LYE3PqJgFC9TAVZVoKxeTrk82nrwxCiHTgMuBJ13OBMwP8i64izwBXTUQDJ1IgadrK9hwjPNwYvMYEYFZaLIeTemgQA3RxIuvQOR1xNAVpbac/aosaSS5TGREU5XQxvA5zsh+hxtce5q+Bu4Fo1/MEnJtyDo+71QKzg9w2nxzdX0lzTRuaTkPTxPEgKKUE6fxTSokQAiHEqGHYpmqft0E7hcMhCQ83MDAQAqPQZ8ZSI5pHH5OSwX9WYrNOTeBq2dtAfl4aBRGjsxOpHqaiTEdqSBZ8CJhCiOHdrPcJITYMH3ZT1G1/TQhxK3ArQGZm5jib6VlyupnDO0p45qevYLO4v2/mSUJqHLn5aR7zFziDqyvx6ogfh0XMiqLdj0k1wWQw6oj/eBpNhkF2nRwsXaISIunpcL9OdKJ1tfWS+noHiVdF0Ko/8aVCpwKmokxL0yyX7ITwpYe5Hvi4EOJSIAznPcxfA3FCCL2rl5kO1Ls72bXH2O8BVq9eHfTZJ3FJMVzz1YuYtySD+659xK8eVVtDJ20N4088ntjehxYeNqGTajzJPD+DN0ye700apEZv29RuNdZQ2kLOn2MwfN5MgyvRgVpWoijTj1pW4uQ1YEop7wXuBXD1MP9bSnm9EOIfwDU4Z8reCGyewHZ6teLchdz6o0/xu++9gH0C086NMmcWjsZTJ7YEQ3iEkVmZccc7t0iJ0GkIDYwJJvZ62WEkVuoJizJNWQ9zWEdTN7ld6TQkugKmyiWrKNOSGpINbB3mt4HnhRA/BvYDTwWnSeN3xc0baKxq5aXH357waxkjDLS2nrp7SLDEb0rh7XD3Q62+aNVZqb86gdzdMZTvqgpiy/xXH3ViIpLqYSrK9KNS4zn5FTCllNuAba6fy/Fhw83JdsO9H6eno48tz38wodex9FvJyDBTV9eJXq850+wN7zk54k+Ec6LL8R1DBDjsJ34WQmC3O3DYHaOGdh1INAmOAH5Hy4z9NJ2l51xDNsfeqwjg3QZGuG78ClTiAkVRpq/TLtOPKdzIXb+9kcbqVg7vLPV+QgBiBy3UVDYTjHmyGqNnUnX8rotFt82jIKw3oHp7dDbeXGtnU3wuZa8c9Vo+PMpE1oJoivePfwbxSJpOo8OVxEAFS0WZvtSkn9M4+XpMfNSEX6O1uoXwKFPQ6hMjHnarHd0fKtjYGkeUI7CsQlZN8n853aR8YRFhEZ7XjuYsTeSvr/ybh377DDfeEe2xnD+yr8yh0TXhx6gLzexIiqKMTa3DdDrtepjDzrhkGTte2z+h12isaCH/3EUU7KuckPotg1aa/nyEebcv4GB4gPdLhWBrYidnfX4ejj+VYhk6dQnO/IUQFu6cJPSpzz7HGWflgZD8e/OZbPt3P1/8uiu70fDIsZCEhVnQNNcQtENDIhDCgZQaUoJlfgmXmgAEGRHjvyerKMrUUpN+TuOAecFnzqCyuI5/PrplQq9jt1i9FxonIQRZV+fyWnjwZuLujOlixS3zMb/UQHvDiXrn5ydw7Y1vjLg2ZM45AsCXvnaUL9xqwmgKNB1gaGRGUhTFTyHa45tsp23ABEjNSprwa/ibLMEfUkqqU4K/RGZ/RDfrrkwjbasJR7+NSz+j56orn0Gn97yGNfBgqSjKdCVR9zDhNA+Y4ZHBu7/oSa9liLRr8hACpF1i7xyi6d3KoNQdGRvOXlM3E7HHx+7oLi5dPYvqzcdYuaZ9zGAZPOo/nKJMV6qHeZoHzJXnLZrQ+oVOwNwYtiWfmMm6IicJghQwNU2Q6DBithvRSUG3zkadfjAodZ/3GlSUl3L1dTrmZO4LSp1eidP6101RTlvDk35mutP2E0xKyeP3PB/0emffsIzOnEgAmh1D/Kdv9GQcKYKXJq+no5+kh08kHTB/LXdcATNc6jhvtwkhoL6gCfPsWOrLnbuYVB4DwdjDrXaZwl8fW0XlEQd6I3z+zj4yM7f73Q5FUaYvFTBP44DZXNPG9s3B7TlJ4NgCA9U9ntco9k3grnHjupspJRcci6Hsw7Ljh+pLT2z59dFuOzZHEgat2u3pA0O5/OD2LA7vOtGL3vWW4KUDmRh07s/xzOBneUVRQkGoZvoRQjwNDG8Qku869gPgv4DhD7rvSClfd3PuJuA3gA54Ukr5oLfrnbYBs+jD8qDXGZ0dT2lv95hlEm0G2gK8zvAWZCdvhdWi93/iTaTUU7N17NR4Q0Nx6MNbcMho+nrTqKlMpK7CyIEPYNvmXmB0r9Zuh+vWz2XxavdD3poOdDoBAjTNOeNWSoiKi+HO3/r9FhRFCQEhOunnT8AjwJ9POv4rKeUvPJ0khNABjwIX4tye8kMhxCtSyqKxLnbaBszW+o6g19nf0M3y6DwO9LgPifGmMCItkrD0WPR65yJ9TYBjOF2eAyQSTRPOXz1XQBzeNUwIsHd0UXGwEoDsZVnIaGcCgciUaP6l9z+n0IbaWMoHT13/aDDomDPHjMmg8cJTBjY/FY/1eDy2uh7DrTtVX7dkz7sDfrUlbe7EJ5NQFGUCyNAckpVSbhdCZI3j1LVAmSvFK0KI54ErgZkZMPu6/Psw94Vj0E7MB22Q7/71bssQlqI+ancUBOV6w4ETYMHHlwLhvp8sJev7zDRtrTnlpbBwA7PCNUp3FANQsC2wdvpq+EuEoijTyzSc9PNVIcQNwF7gm1LKk3tQs4GRH461wDpvlZ62AXP1Bfn8/TdveC/op4bXSph35lqO9Zy6vVZWZAy1fwz+NQHsuWbSLSfuYmoI9FIgJEhNkD4YxpBOktihYa/opfVYO+097afUM29eIrUHyqnq9723qtNrZMxPITImDHD1OcWJPbU9O953BiRxSTE+X1NRlNAyRQEzUQixd8Tz37v2WB7L48CPcH74/Aj4JfDFk8q4ezNeZ2yetgFz4epsElLjAtog2h0BpPfoOObmtWO9XZz7gwup++lWHG5Sz41X5pI51DxTSLR0/++5eH0OhcXOL0ue7rBKKVmyKIWmY40M+hEsc1bMobmsnmO7j/jb7FPMWZQecB2Koky+KZz00yqlXO3PCVLKpuGfhRB/AF51U6wWyBjxPB2o91b3aZscUKfXcdG1Z01I3Y4DLR5f+09EL/ZfnU9klplg/X5FpyQcnwjkjqZ5/2dMTo6heMcRmmt8m5IUlxTNwhWZFG87TFttoNOYnOy2yUiOoCjKRJBSTPpjPIQQqSOefgJwd4/sQ2CBECJbCGEEPgu84q3u07aHCbDmwnyee/iU2cQBCUuO4siZUTDQ77FMRW838d9YxpDDzlxTFNFDYKrrp/oX28Z1TZ2XkYKuZu+5Zs3mCJqsvgWsnOVzKN9bQkFZnU/lfWXz8fqKooSeUJwlK4R4DtiAc+i2Frgf2CCEWI5ziLUS+JKrbBrO5SOXSiltQoivAm/iXFbytJSy0Nv1TuuAmbcqm1XnLWLf1jEnPvksflUaTVel0eJlaQlAx5BzKUaB1XWveRacc9Maav/0od/XdTjGXoHZWNWKIS4GqyWwgJSQEkdCYgTH9hzFMhCMXT5Hs6uAqSjTkgzdWbLXujn8lIey9cClI56/DvjVozpth2TBuY7xi/ddHbT6ItalUelDsHQn1mSi5bXicZ3b3dRJ/CzPE2Ysg1bS0uK81DJ2LzUrL42Blg6KdxyZkGAJakhWUZTp7bQOmABz89NZf/mKoNRlixt/ppquoSFMs8a3DrG6oIqU5Gh0es//XNFRgW2dNdjZQ097gHtueqGGZBVl+pou9zAn0mkfMAE++/VLglJP+2/3kRMdP+7zY65dhtCP75fg8LsHsTa1kJAS6/b1mlr39zGNJj2L8tMIDx8dUCNjT6zpzF83l7qjwb1fqWmCrCWZxKaYMUVHYIqOQGdS+2EqyvTknCU72Y9Qc1rfwxzW7sOkGF9Yeyyk7umke1UkrYMD2LzcWzzZdmM3Gb+5EPNv9tNT4nmmrScOhyTRHE53ex/Wk/bhzEiLoavTORHJaNITHm5gdno8R9/ay6HDJQCYU+OZvTCd/s4+msqbyF87H4Tg8DuHTrlWREw4qfNSMJgMlB+oxDLofph2yTmLaahswRhmQG/QIaVkoG+I9oZOqo82jirbzyBWiw2DcUb82inKaSUUe3yT7bT/5Co7VM391z0atPoaXjlKzCuQ/sOz2NvnOQm7J7V9PaQmRY0rYALowsNOCZYA1oEh0mbHEhtpoHRfBV1WO6KhiaHeE3lg2xs6aG84kfDi4JaDo+qIT47FOmQlNjkOITQqCp29zlmZCUTHhlO6b3R+3uylcyjZX4HVxzWn1iErd5x9P9968kvMXZLp83tWFGVqTcNMPxPitA+Y0fGR6PQadtu49vrwKHJfO8bFOix2/+7LSaDn8wtZvDSDvverhtPJohtOG+d6LoQzlc7w+kspJVJCb4/7BOydLT30dvZR3z1AdHwkyelmGgorfWpTWKSJeSvm0t87SFVRHf0Vo78INFe30SIE2cuyiYoNp+ZIHZmL0jl2qMbnYDms/FA1X1t/Hzc/cC1X3n4hOt2MuCugKNObdM6UnelO+4CZnJHARdet599/fi+o9TZsPsJZEfn0Loz3q6ep1zQi/17F0SPjSwaQmGbCaDJgGbKOOt5UfaINKRlmit89ePKpbi1en8eRvRUU73GXu+gEKSVVxc4epxCCgp2lfrb8BJvVzu++9Vc+PFDL1+67irTMhHHXpSjK5AjFdZiTbUZ8vf/qz6/lhns+HlAdmk7DnBzDwlVZ5K2cQ+6KOUQVdZPyz1quPBbGsmjfPvRtDgcxV+cQHR9JXKL/s2Y7W3uYt2wOeWvmeSxz8rZgYxGaNmYWIXf8Le9Oztr5FB+s5sufepTX/rEnKHUqijIxJGqWLMyAHiY4U8d99q5LaGvq5LU/bh+zbNLseIwmPULTnNtwuYLPUGcXtUUVNBdVnHrSu4AmOPeLa/koP4wey9jrGD+0tzNHOrC0DbJozVy/9u60WewUf1SJ0AQx5ki62/t8PvdkS85dROEHZd4LBllyVhL1gzA4ZAfs/PZHr7Br6xG+/j+fICEpetLboyiKN6E5a3WyzYiACc5e11ceupashbN5/J7n0Rl0xJqjCI8KIzo27Hi5+qIqGmr8n8yDQ1Lz5G5Wff1stqWNXXTJ+/20djkn4xTtrSAtK5H6Sv+uKR0SO4LwqDCyF82mo7mbhkrfJxIZw4y01gU3Mb0vktITiMhOpfmkZTAfvl/CbVf/L1/7/v6FTSkAACAASURBVJWcc5GH/dMURZkyahBoBgVMcAbNy79wLovXZvPbrzzJ4e3FCCGCOhxY9ev3WPXIZezrP3VrLYBzdIk07tw76lhsYrTfAROgr3uQhORYij6qYvXZC4iJjwAJkVEm5izJZLB3iP6ufvo6e3E4Rr/HucuzKPmo0u9rBiL/vHxKm/pp9bBmtKdrgAf++3l2X76c2++5jKgYP/b/VBRlQoXiEOlkm1EBc1j24kx+8M9v8fR3n+PNP74b1Aw0ArD8fCcbPrOcQ+nQPjiIUadjXmQss3oEjb/ed8o5R/dXkZKZQGO1/xOB2pq6mLcwjQPbijzOWBVRUURFmYiIDic8KgxTmAF9mIHw6DAGegbdnhMsUQnRZJ+dj2VwiMJjvr2/d149QG1pPTfdvpEVG1VvU1GmmpQqYMIMDZgAMQnRfP2JW7n4C+fxjbO/H9Q8p31VHfT9bCtn3XI2nf3Q+kEl9iE7DR7KOxwSnV7HrHQzzbXue6ZjGRqyEh4VhnWo12OZgd4hBnpHL0nJP2M+h7d5TdAfkNRF6RQU1hOXEOnzOVJKLNWN3H3hj/jEHZdwy0+vwximsgQpylRS9zBnyCzZsSxct4D7/vHNCalb12CheVsljiHvwbiuooXYhCh0Bp3f16ktb6HfBvOW+pcMoLasCaFNzH+CxZesZNGVZ9A8jg5smElP2V7nMpeX//ff3L76Hkr2+T4xSlGU4JNy8h+hZsYHTICzrlzDHY/9V9DrHej3b9eP0kM1LFyVNa5r2a122jsHyT87z+dzOlt7WLDa8/KU8Vp04XJKKjsoLmqgu9PzvqGeDFns5J+35Pjz6uI67jjrezz/s81etzpTFGViqGUlM3hI9mRX3HYRfV39PHXvs0Grs6WonMWrcyncX+PzOd3t7gNMRHQY8Ukx1JU3H38O0D/iHmRnaw/pmf4lh7e7iT855y5GFxuNTgMhJa0ltTQUe07OHp+eQOrCDIROw2EyUlzkafDZd01DEBkbQV+X8+/DbrPz1L1/Y/drH/Gtp28nbV5KwNdQFMU3ktAMYJNNBcwRPv2tj9Pd2s1bz2yjqzXwra66mro4+Noell60kqIC34JITVkTN9x9Od3tvbz6zHt8+qsXcvWtGzGFG9EbdHS192IZtBJrjqKiuI6vX/7L4+fmLs2g/FC1X20sL6hh5SfOYEgKBi02+vsslNZ3Qv3I+6F68j5+BvT2Yxu0EBYfidDpGLRDS0sPnR39dFZ2jKrXbI5gaMiGTqeRmhxNU6t/60XbuwZZvH4Rh14fPaO44P0j3Lby23z5Vzdx8Rc2+JWkQVGU8QvBEdJJpwLmCJqmcevPb2DTzefzy5sfo+iDkqDUW/pBMXPXLaK8qN5r2diEKD5520aMJgM3fvtywiJMo183n8gONNzLHNbT2UfW4nQKP/AvbV175wBVtWOvyTxaPCLgtwwc/zF9dhxps6KwWOwYjTr0UtJe107Dkdrjy3V6SiEuL92vNgEMehh9Hegd5Jf/9QS7Xt3HN35/K7GJnjfXVhRFCRYVMN3IzJvNb3b8hIbyJh78/P8GHDgHegYofecjMs7Ip75q7PWWK87OxWhyblR9crA8mRCjb0HXV7djHk+mHD97aSaTnrlzzPQ0tNNSVE3doNXrOWEDgyxbmMyho81Ix6nfVYUATWjYR9yjNA65TzQ/bMfmDyneXcpdf7iNdZcGZ5NwRVHcUMtKADXpZ0ypc5P5/j++yeL1uQHXJaUkyuTMSTuW917dzyFfE5u7m0Y2jnETXxM3CGBJThJaawdH3iumrqwJyxjBMizCSP6abHIWp9FUXMWBzbuYZellfmbc8TJps6LITTRirK0nutO5TtNo0LEoNYLGUu+bWrc3dvK9Kx7kf7/yJIP9YwdYRVECIKfgEWJUwPQiMc3MvX+9k7ikwIf9irYdJnfJ7DHL2Kx2HvraM+NKphBrjqSjaeLS3c2ZY6bgvSMM+jD7V2/UERdl4NDbBzmy8wh21/tpLG+m9O39zI/VyI6Q1L1/mCM7jzLYN0R/dz/zMuNIToigZGcxrX6sSf2/J7bw5TX3cOxg5XjfnqIoY1CzZNWQrE+S5yQRbY6is6X7+LH5G5YgNIGUULb1kM91HX59D3kXraasYHTvKX3eLPJWZiElaJqgvbmLWbPNY9ZlijBy5qalzryyQxb2v1tAp5/7U3oTHm4gKyMee+8AFQd8Wwup02vMzUnmyI4jHsuUfnjqdmJD/RbK3t5P/lk5VPq5JAeg5kg9X133Hb74k2v55DcuQ9PU90FFCZZQXBc52VTA9NG1D93In59y7nRitzs41trL8JjB4s+dT+Xre+hr921mrclw6gf5wlXZ3PXw9X61KSktnvuecq4fveP8H/u8mbPeoCNrYRomk56miiaiDScVkJLcnGQ0i4XW2naOvFfsc5t0eo35eSkUjxEsvXG4W+viI5vVzu/v/it7/r2fbz55GylZs8Zdl6IoTsPbe810KmD66MKPr0RnMvDzH27GftIHelFJM3FrFpFmHaD0Xe+9zZLth8les5CKIyNmngb4uzg8UUgIQd6KTAZ7B7DYwGqxOXdkiQnD4XDgsDlormykZOeJIKhpAmmIRG/QYzLpyU6Npuh934PkSAtXzOHQ275tXu1JeUENUfGR9HaMf+uyA1sL+dLyu/naIzdzwefODqg9ijLjSSAEA6YQ4mngcqBZSpnvOvZz4ArAAhwDviClPOVelRCiEugB7IBNSrna2/VUwPTDxouXMDRo5Vc/ffWU1zo7+ukEFl1/Po72Lqp2FjPQ1YfQaciTAuxAzwDhhtG/fGEB5EqtLK6jorAWgMVr53LoHf8CVmN5E/kXLKf0cC0WKSmqGF/igYUrswIOlgApWUlUF9UGXE9/zwAP3fgIe/69nzsevZmoON/z2SqKMlqIDsn+CXgE+POIY1uAe6WUNiHEQ8C9wLc9nH+elNLnraJUwPTTRZcv59k/vkdzo/stqopLnZl4tEXzSE6IprOrn/nmMAo3f3C8TO66BcTPTmCe1JwzVCWkZieOqz0/+cITvPcv5+J+vUFH3RHfswqNZO3uDWjXlsycFCo+Cs5m1JWFtUTEhhMdE0FbfTsOh8RmGf+92a3P76DogxLu/evXWHxW4DOeFWVGCsGAKaXcLoTIOunYWyOe7gKuCdb1vAZMIUQYsB0wucq/KKW8XwhxPvBznDNte4GbpJTB+cQMYTqdRrw50mPAHOawS5qbnZOEihp6WXrtBg49tw2AyJhwvvuHWwLOUlOyv/J4sDSnxKKzW2mqaB5XXaV7jxE/P53OVs87nngSFmGkq7aZ/i7/88Z60t81wGDvEEnpiSSmxVG8qzSgoNlU1cJdG37Atfdcxee+/0n0BvVdUVF8F5qzVn3wReDvHl6TwFtCCAn8Tkr5e2+V+TKNcAjYKKVcBiwHNgkhzgAeB66XUi4H/gZ8z5fWT3dFh2s46kPGnpMdPtZGwiVnkP+58xmKiODl323xe+PqoQELNqsNq8XGB6/v528/f5Ul63PJX59L6rwUulvGn87P4ZCkpcX6fZ6m00ieFUVHQ/CXszjsDpqqWijYWULa/BQS0s0sPXcRpojxDV877A6e/clL3Pmx71NX1hjk1irKaW5q1mEmCiH2jnjc6mtzhRDfBWyApwTh66WUK4FLgK8IIc7xVqfXr9nS+ak+3O0wuB7Db2d4cWIs4H8UmYb+8dcPvBfyoL2tj/a2PowmPUOvHOTZX77BnNxUABJS4rjmKxeSnJHAuy/uYd+2Ino7+xnoG+Siz55FYlocb/99F221bfR19dNc287cJRk0N3TS1zNExtxEBvsC2wy6ubKJxefnU/hOgc/nLFyeweF3Dwd0XW+EENSWOgNcR1M3C1Zmc3S3f+n/RirZW87tq77NnY/dwvnXqwlBiuLV1GX6afVlMs7JhBA34pwMdL700DORUta7/mwWQrwMrMU5muqRT+NSQggdsA+YDzwqpdwthLgFeF0IMQB0A2f4+mamqx3/OcKO/4x/uQQ4Z6TOTYykdMdRAIr3Vhx/7f1X97s956kfvXz859SsRMyzzejDTVSUNB0/3lzfhTk1jrb6DndV+EYIDlW0s3hjPsd2HBlzmUpkTDjZuSlBmeTjr7IDleiN+oCGaAd6B3nwhkf48I0DfPW3X1QTghTFmxC8h+mOEGITzkk+50op3d4nEkJEApqUssf180XAD73V7dPKbiml3TX0mg6sFULkA98ALpVSpgN/BB720LBbh7vTLS0tvlwuZEVGnZrb9YJLlxIWfvJCRs/yc5OPB8vxaKhs5eiBahqq20YdN4UZAt5hZTjHa2FlOylrFhAZG+62nKbTSE2LnZJgCTB7fkpAwXKkd/72Pl9acTeH/Vhrqigzk5iCh5cWCfEc8AGQK4SoFULcjHPWbDSwRQhxQAjxhKtsmhDiddepycD7QoiDwB7gNSnlG16v5+99NCHE/UA/cJuUcp7rWCbwhpRy0Vjnrl69Wu7du3esIiGvoqyJ0iMNWCw25sxNYsnyORzYV8HdX/mLT+evXJbBwc0fBr1dUkrmZZsp3XdqBh1fJWcl0ZR4YqF/SkIkjtoW2k66PxkRHYbRYaO93vfUdcEkBCw+M4ejH5YxNI6MQO5omuD6732S6797NTq9Lih1KtOHEGLfeIb+ZgpTdrpM/cHXJv26VTfdE1L/Ll57mEKIJCFEnOvncOACoBiIFULkuIpd6Dp22suen8xFly/n8qtXs2T5HACWr8rmvp9+iow53peGWDXBHb+4joylmSy8eBkLr1hF5oqsgNulaYLWWp+XE7l18nenxrY+mJ1IUvqJFH2ZC1KICddNeLCMTohi7hfXkbEs0207C3aWEJcSR/7H8oJyPYdD8pcfvsg3N/4PTVXTeyREUSaESr7u0z3MVOAZ131MDXhBSvmqEOK/gH8KIRxAB87puzPWx85byIK8VJ75wzYO7qtkoN9Cb8+pk3AG+q1c8rmP8VF5K1vfLjx+fPGGRZRsKxr39U3hRtKyzHQ0jb3cZSwO+6nrMJvb+4mOj2FOTDiRkUZKdh7BMsG7gqTmpmI4L5t3EyxEXJrMurWzGQjT0L1TSV3RiRy8LTXtJGeMb/2qJ4U7jvKlFXfzjSdu5dxPnxnUuhVl2grRTD+TzZdZsoeAUzYblFK+DLx86hkzV3JqHHffdxUAfb2DvL75I44UOj/g585PZvGyDJa6epObLl02KmAWNXSx+Px8qnaXMdjr32xXCcTONhMWH+W17FgsA+636urpszCo18gJMzLUP+R3Fr+sdXMxpsfSsrPC6/KTsCgT+65Ox6I5h1r7NcnWWc52RV6ZxsKLsuh8dNfxJAt22/iTLXjS19XPj6/9NXvfOsiXf30T4VFh3k9SFOW0p1ZvT5DIqDA+df1ZHl9fuSqLtevmsWe3856jlFBQ0wFpCRiNemJiwkkx6Ti6tfCUcyWQsiAF89xkrAYdTa291Lb3YR4IrOdns3qeSGO1OSis62HhZ86h6o0PGegaGLMuvUHH7JvXoh9ysCvFQZ8m0eXM48wuEz0vFtDV3E1kbAR2ux3bkI3sjXl058ZxKNKCRXM/FtOnSfZGDpH39XVEPH0Im8XudjPqYHnjj1s5/P4RvvPsHeSsmjth11GU6SBEU+NNKhUwp4gQgjvu2sQtN/6ewZM2YbZYbLS29hCRYSYmOZaelm5S89KIy0xkSKejvqmL2u5BaisCu2eZuz4PY1Q4DqMRG6DvH6Cwa+weW3FJM0ln5GNuaKbuUNWo12ISo0mYPwvbymQ+iLNQJka/L7smeD/eArcswCz1NOHAgIYVB2WaxJkjY2wLh4yY/nyYxIwEjGFGyg9WsmDlXAb7h6g9Wu93Mghv6kobuHP999SWYYqiAqYKmBOtf8DCsy/swmZzuP0wX3tODtvfOrUXCVBd044xKQ5dipmqAQtV5YEFSIDI2AiyV8+nyqpxtN8KQxacSf1919LcgyEsisUfX0fBK7uPD9EaPr+UrRFD3usTgnbhDMwD+D6kek6LnrqnPqDfYqOjoRMhBInpZsoOONeyLl6fw5HdZQHlxHVneMuwfVsOcfefvoI5JS6o9SvKtKDuYaqAOdGMBh3PvrBrzOGMpSsyKd5f7fY1iz/rDb3kpp2/Zj6NYZEUdNnAj0DljtVqp6Cmi0WfOZeuw2XYLsgicncTbIj12g5/GRywoWCIkn/tHnVcSklLzYn1qIU7S1jysTwObZ+YCdv7thziS8u/xbf++GXWXnLKbX1FOa0J1cP0LXGBMn79AxZ0urH/mosqmpizIAgbHY8RlZdcuooyh5HegeAs+B9W29TNuf91ESVmydbz4oIeLGfZdZzxWiMl/zrgU/mq4joSRyyDCbbOlm6+e/mDPH7XM1iDlDxBUULeVCwpCcEArXqYE2z7jhJsNseYZWw2By19g8QnRNHR5v9uIWOJiA4nfcNSDjf0Bi2YCQEr1sxl0xXLOeucXIxGPXJXPA9uGTMNo9/yB41ozxyi2o+lMt1tPazYmE9r7cSuE33pN687JwT99Wuk56RN6LUUZeoJNSSLCpgTLibafXq5k3V1D5CZbkbfM4DNMr7hUnFShpq0+SnYsmZT0hCcIJyYFM3Fly/n4iuWk5I6+j7e59Ys56WDRZQ0B36fFSDFrqf/4Z1+LRvR6XUsPDOH/VsLmZOfQWRMBO0NHTSOc8szb0r3lXP7qnu447FbuPDzXjc6UJTpLQR7fJNNBcwJljUnweey1bXt5C+ZzZF91X6vdQQo7LGRdcU6Wt4rICM/k3K7HktnYDuYaDrBuvULuOSKFaw5Yz46vfvhZYNOx13nree2v28O6HrDIhzCr2AZmxiNOc1M4c4ShKZRc7QBgJiEKOYum0P5wSovNYzPYP8QP7vpUQ5uK+T2h28kMjZiQq6jKFNOBUwVMCdaWkocOfOTKSlr8l5YSiqrWjljVRY1Ne2YDDr6BixYrHbiYsKJRmC32o/nJBZCILQRodX14+JPn82efdUgxh4KHkt6pplNV6zggkuWYk7wLSHChgXZXJ6fx6sFge3oMsuuQ/er3T5PS8rIS6Ove4DKwtpTXutu6z3e8zQY9bQ3jN7NRQiBwyGPz2DWNOdzIYTbWc1CCDRNHP+71uk0dAY95Yer+eUtT/C1R75IfLKaRauchlTAVAFzoun1On7z0LU89uRW/u/fnnf30GuCBZqRxrJWPioYvbmxpgn6/FigPycnZVz3K01hBjZelM9Fly1j0ZJ0hJ91CCG4Ij834ICZ/LvDdA+6zzp0skVn5VC6v2rM3Us6mrroaOrCYNTz0BvfYfFZuQG1T1FmHJUaD1ABc1JEhBu56bqz2Lr9CL197hfnx0eGUbX/1B4SOBOD+6O+spXFq+ZQXt3GwKBvMzk/ee0ZfPaG9cTGBTakuC4rgyxzHJXtY6fAG+ncFj2Nf9l3POh1+zD7VKfXsejMBRzeUeJzYF932QoVLBVlnNSyErWsZNIkJkRz3afWeXw9ycfJQb6wWmyUfViJ3e77b7jDIQMOlgDhBgOPfeZKwvTev4vpHXDhoSEqH9/BYO8gNovNp30ujWEG5izOoGBnqV+94MrCWrUURFHGSy0rUQFzMpnjIz2+drShneSM4K0ftNnspCbH+FTWnBjFueePuZWpX+Ylmrn7grPHLJPo0HHm600+r68cKWf1PCoKavw+r62+gz/d/wKv/uHtoKfQUxTl9KeGZCdRQXGdx9fsEmRmHMbmbixDgfeC5q/N5mhZs9d7mXmLZ3PfA9eQOMu34Oqr69cs57l9hyhtaTvltYVDRkzPHKa60fdh25HEOPO5DvQO8uKvXgMgLimWj121Zlz1KMpMpIZkVQ9zUjU1d4/5erfFxqObv058UnTA1+obsHoNlpd8fAW/ePSGoAfLYQ9ccREr0lNHHTu7TY/117voGGewdAr8f25dWaP3QoqinCDF5D9CjOphTqLVK7L48KPKU46fuz6HizYuxmTSkz43iSde/QbPP/4uL//p/XFfq62kkeiUOHrc7K2p02l8+a6LueLq1eOu3xdLZ6fw/Bc+S2VbB195fjNzd7Rw9B+7vZ/ohTUIPfAgZ/BTlNNbiN5TnGwqYE6i887O40/P7mDAtWTi7LMWcN+3r8BoGP3PEBMfya3fuYLe7gG2vLRvXNeyDNnISo3B2m0kMi4Cu9VOe7+F/iEb9/30GvKXZQb8ftwZ6Buk4VgTVa7h576ufurLGkh9s4CjRe5nAfsiOj6SaHMUibPNHHr/qN9LXk720TsFfPqbVwRUh6IoM4sKmJMoeVYMv37ws3zzey/Q2zvE/kPV9PQMkmB2nxjgrgc/zYYrlvPvv+9h1ztFfm9bVbKrfNTzsAgjP3/+dubnBTf3aVNVCzUl9Xz47wO8/uQ7DA34t12YJ0vOWchA7xA6nUZTbSv1FS00VLYGHCwB9r9bwL2X/ZSENDPX3XMVKdlJaq9LRRmL6mGqgDnZ8nJS+fLN5/Gz37zB0KDN604mK9fnsHJ9Dr3dAxTtq+SlP77HwV3HxnXtmPhIMucGYVcUwG6z84+HX2Xn5g8p2VeOwz7+rEKeDPYNUXagEk2nIV3Zd4Lpo3cKANjyl+0sP28xX/rZ55i7ZGJ63ooy3alJPypgTolNF+Sj02lYbXbifMw9GhUTztrzFrJ8/QL+8+pBnn1kC021Hd5PHKG1sYvy4nryAhyO7Wrt5oHP/Zb97xwOqB53NJ1GbGI0sbNiaa3rcKao8zNxw3gc2FrIV8/8HpfefB5f/c0XJvx6ijLtqICpZslOBZ1OY9MF+VyxaZnf5xqNei68ehX/87ub0Bt03k8YYe2GvICDZUVBNXed/0MG3EwmClS0OYpFZ+YgNUFlYS2dLWPPKg42u83Oq79/h+aaU5fCKMqMF4KJC4QQTwshmoUQBSOOmYUQW4QQpa4/4z2ce6OrTKkQ4kZf/gpUwJym0uYk8v1Hb2BW2uhE37HmSBYtTCZ/aRrhkSbCIoyERxqJiDKhaYLe7v6Arlv6UQW1Rxso+aiClOzgDO8CxKfEER4dQeEHpXQ19wR9+NVXUkp2vvLhlFxbUUKVkFPz8MGfgE0nHbsHeEdKuQB4x/V89PsRwgzcD6wD1gL3ewqsI6kh2WnKYNSzdkMe81/4Cn//3Va2v36I3u4BZkXrKXjT/czaj2rrGbj3MqJi/EuBZ7c7+NP9L/D+y3vo7xkAQDokUfGRiEr3u3r4KzV7FsV7xndvNtie/v4LnHvNmcQnx051UxQldITgukgp5XYhRNZJh68ENrh+fgbYBnz7pDIXA1uklO0AQogtOAPvc2NdT/Uwp6G+ngEevutZ7rj05/zqm8+i9fay8YJcMhLCOPpeocfz+rsH+M4VD3HskO97Q3a39/L9K3/GC7/4P+qPNdE5IvlC/bFmwiJNAb0XgNy18ynaXRZwPcEy1D/E5sfenOpmKEpomZoh2UQhxN4Rj1t9aGmylLIBwPWnu6Gw2cDI/Jq1rmNjUj3MacThcNBQ2crXL/8FvV0Dx48P5Kdy+D+eA+VI1cV1PHD9b3nkgx8THhU2ZtnKolp+cM3DNJS738tzoHeQ3NVzObK71Pc3AegNOpIyEoh09XQ7mrqnbAjWk47mrqlugqKElCmaJdsqpZyIDCvuPnC8vkMVMENURXE9u7ccpq2xi+a6Dlrq2mmoamWw/9Q1jgWH68n92CKq9h9j0MP2YSPVljZw3yd+wVVfvZgzLluJTn/q5KGdr+zlZ1983OvknpJ95Sw6M5eaI7X0dPSd8np8ShwpWUloOg1N0+jp7GOgd5DGqjYgNCfXGMMMfOKrJ98WUZQZbvrMkm0SQqRKKRuEEKlAs5sytZwYtgVIxzl0OyYVMEPUYN8Qzzz0qk9lhRCUFjcyOyeDvoZW2n3I03rovWIOvVdM4ux4Lr5xA5fevJHE2WYcDgd/e+Bl/vLjl3y6tpRQvKeM3NVzMZj0NFY4fzdTs2cxOGClbH8FnS09PtUVKiyDVmpKGshanDHVTVGU0OD7JJxQ8ApwI/Cg68/Nbsq8CTwwYqLPRcC93ipWATNE5ayYQ4w5iu72Xp/PqatuJzIqnLxzZnP0vSKfJuO01nXw7AMv89xDm1l/1RqMJgPv/M3/HLZH947OKtRaH0hy9cCl5s4mbnYCSImUuP6Ux9d0yuHjgN5kwBBuHHX+oY+qWbNpOWEnHVeUGSsEA6YQ4jmcPcVEIUQtzpmvDwIvCCFuBqqBT7nKrgZuk1LeIqVsF0L8CBieEv/D4QlAY1EBM0TpdBprzlvIO//0b4lDX+8QOk2g02t+pdJz2B2898/dJKSZmT0/Zdrv5pGUk05BQcO4zy+v7mL9ZStYfsa8ILZKUaaxEAyYUsprPbx0vpuye4FbRjx/Gnjan+upgBnC1py/2O+ACVB6tIkFZ+bR09RBa22bT/c1h7XVt2MMM7Bw3XyKQ2jm6skyl2ejNxiw26w4bA4cdonDZgNNw2G34xCBTQCXUqLXq0nkijJsGg3JThj1iRDCyg7XeC/khs1i52hhA/Wtg+jM8SSkeV2PO4pl0Erx7jKWfCwPTQut2avD+g3hVLYOUtNpp65X0jAATVY9TUMaLTYDRYc9b9bti8F+C48/8H/0+/FlQ1GU05sKmCGq5GA1Lz7+TsD1DPRbSMv1urzIrcPvHyFrcQbRHnZTmUyaJkYtPZmML7vlxQ08/uNXJuFKijINhGBqvMmmAmaImpUeT/ysmCDVNv5eYvnhagxGPek5qUFqi39MUWHkXbYO5mQgsjLQsjPRz8+io9X3yVCBaKxpp6rU/TpURZkxQjc13qRSATNExSVEk7diTlDqOlbWQlZ+BomzzeM6v72xk6aqFhauWxCU9vgqKTuZxBW5lBSNnrwzEVuJeVKwr5KHv/vipF1PUZTQpQJmCDvz4qVBqWew30J1Qx/9wjDuhOnWnjXXSwAAIABJREFUIRvFu0snLWiuvWQ59/z5q8SEwHBwZ1uvupepKGpIVgXMUHbmpqVExoQHrb7BASsDwhDQ8KrOMHG/MmGRJpLnJHHLA9fyw5f+m/w1c4mKGTt932Roru/kV999kWPF9SpwKjOXCphqWUkoi4qN4JLrzwrK5J9hPV0DDBr1LFqfR9GOI0Gr1x9X33kJYeEmNN2J4LtgZTYrz8/HYDKMmtxzziVLMYUZ2fm2b7lyJ8r7bxaw461C7nvkc5yxcdGUtkVRJpsgNO8pTjYVMENc+vzkoNdptdgoLm4if+MSiv5TOKn3BAFu+P41XhO/D7vo6tWER5imPGCCc23mGy/uZcnauUT62H5FOW2ogKmGZEPd7LnB26R5JCEERQfrWHBm3vFjmiZOScSuNzq/U8UkRJG7eh49bYHPTp3sAB1Mu7cWs3d7yVQ3Q1Eml5olC6geZsjLmJ9MVFwEvZ39E1J/S3MP4VFhzFsxF2HQ0drSR2JSJEd2FGMdspE2L5nYFDO2/oGgZf6x+xEwLUPWoFwzWMxJ0Zy9KX+qm6Eoky8EA9hkUz3MEBdrjuKRN+4mN0hLTE7W2daHFh9HUVEjhQfraKrvpPBgHdlrcog2R1FdXIdDSvThgW8UPcyfvS/vu+0ZHvjG34J27UC1t/RQXxWa25IpyoRSk35UwJwOkjMSuPfxL0xY0Bxws8dmWXEjuvg4Fqyeh6YJCg83sGh97oRcfyyhmJjvjRc/xOGYvsPKijIeakhWDclOG8kZCTz8yl2UHqzm9b/u4O0XduNwTOxvVHdnP92d/UAHqenx1BQHlp8VYP1VawiPCl5vdSr88+n3SEiK4RM3fWyqm6IokycEA9hk89rDFEKECSH2CCEOCiEKhRD/4zouhBA/EUKUCCGKhRB3THxzZzZN08hdkcU3fnk9P3r2yxjDDJN3cQF9nX3jPl3Tadz60PV8/7k70Rt8/57W0jC1+2p6kpAcrLSFijINTMVwbAgGaF+GZIeAjVLKZcByYJMQ4gzgJiADyJNSLgSen7BWKqdYeU4ed/7M01ZwwddQ00HOxxaRO45MP+bUOH7+1vf45J2X+nX/EpzrRkPRseL6qW6CokwqNSTrQ8CUTsNrCQyuhwRux7lLtcNVrnnCWqm4tfGTa7jshskbFiwpbKCta4icNb5vqrzs3EU8tusB8sd7/9PPADtZtv/7EFKG4P9oRZkoqofp26QfIYROCHEAaAa2SPn/7d15fFTV+fjxzzOTfd8JEHbCvsmiUFwQV1xwqeKudamV1qqtttXa6rfutdW6Vf3R2sUquLTugjugKC4sEnYIYQ1LgEAWErLN+f0xkxBCSCaTmXtvMs/79bovknvvzHlOwtwn59xzzjXfAP2AS0RkkYjMERFrV+ZWAJx9lbX30Yp3l7N+QzEDjx9CzyE5LZ57+Z3n8/Dsu0jtkhxweQ7Nl+zcto+8bwvsDkMpy2gL08+EaYypM8aMAnKAY0VkGBANHDTGjAX+BvyjudeKyI2+pLpo9+7dwYpb+fQZ0p3Lf3GmpfczRYT1q3dStK+aEScPI3dMX5IyEhuOxyfH8Yc3buea/7sYt7t9A7Gd3Ih7+v/e0lamCh/awmzbtBJjzH5gHnAmsA34n+/Qm0Czj9Ywxswwxow1xozNzMxsR6jqaK6642xeWvwAIycOsLTc6upaVuRtZ8Om/cRkpZGSlUzfET155qsHGH/W6OCU4bCFCxor21/Bkq/W2x2GUqFnR7LsiAlTRDJFJMX3dSxwKrAGeAuY7DvtJEDXC7NRYkoc1/52qm3l79lVRlqfbIacPJKouOBNGzEhnjrTHqX7K9iwSgf/qM5PbNqcxp8WZldgrojkAd/hvYf5HvAI8EMRWQ48DNwQujCVPwaO6sVZFt/TbGxT/m5mv/Ql77+4wLYYrORyCUNG97Y7DKWsoS3M1hcuMMbkAcc0s38/cHYoglKBm37/RXTJSeXlv3xA9UF7ujM3r9vBnh37yeiaYkv5VnFHuBk4vOWBT51JRVUNYIh0u4lsski/UnYQkYHAq4129QXuMcY80eicScDbwEbfrjeMMfcFUp6u9NPJRES6mXbz6fQb1oNXn5xDVWVNw0hTEWnUzyGIyOGjUMW3/7Dv68/FOwKnuWGrwqG/BgXKist4+bH3SeyZ5e1aafSaujoPxhhcLu8+Y8Dj8eCpO/TnpLgOxTZgZM9WB9YYwxFD6gRpNlzD0bt6jDF4DLjkUMzmsBcaMIcqGxkVwcv/bx4xsZEs+24j1dW1h72fy1eP9D5p5KfX4Xa5cLsOle4SaShHGhXTuL4i0qaBRSKHfn8u389QfOXUl9e4/vU/D48xh+rcpLz6/XNX5FNZXUtWUjyPXXMOo/p08zsu1fE5cdSqMWYt3vUBEBE3UIh3TE1TXxhjzmlveZowO6kxkwZTV1PLvVc+a0v5G/K2UpGY2PqJHdzSRZtaPScXD1/vLQt9MBYpKj3Atc++zm/On8QlPxjR5sUoVAflwITZxCnABmPM5lAVoIuvd2LjTh1GSoZNSUunWzQwaVF2hxB0tXUeHvzfZ/xu1occbNKyVp2U8+9hXgrMOsqxCb7lXeeIyNA2v7OPJsxOTES45i77Rs4qr5qt5aRGdL6kCfDOotVc9fQrbNtbYncoKpSMbQsXZNTP4/dtNzYXnohEAVOB15s5vATo5Vve9Wm8MzwCogmzkzvjiolMmDLS8nK1m+6QyPhI9tVU2R1GyKwp3M0lj7/Ml2s22R2KCiV7Wph76ufx+7YZR4luCrDEGLPriLCNKa1f3tUYMxuIFJGMQH4EmjA7ORHh9//8Cc9/fg8/uf9ieg+yZqBGRflBoov3E1W8j/6ZsZaU6VSuCJdz1/gLktLKKqb/7U1mfBz6x84pezh8abzLOEp3rIhki+8veBE5Fm/eC+gp8Joww4CI0GtgV86/cTLPfPpbJpxpTYuzuqqGmqrasF8+rrYiPO7xGQNPz/mK2/71LuUHO2+LOmw59B6miMQBpwFvNNp3k4jc5Pv2ImCFiCwDngIuNQFelDRhhhl3hJvpD02zdDDQ7m3F5PZNt6w8p6mrqbM7BEvNXbGBy56YxYadAf0RrxzKqS1MY0yFMSbdGFPSaN/zxpjnfV8/Y4wZaowZaYwZb4z5KtCfgSbMMJTZPY2nPr6LHgOyLSmvbN8BChauZfjATEcudxVqnbw3tlmbivZx+ZOz+DQv3+5QVDDY0bp0YMeUJswwldktlcfeucOyBds9HsOKeasYOjCTYQOzGDowk+jo8JgG7I4Kj3o2VVFVw23/epenZn+p9zU7A02YunBBOEtMjee+WTez+LNV5C/fwoL3lrJl7Y6Qlrly3qqGr9OG9qKqKgzu74X5Pdy/ffIt63fs4eErziQhJngL8yvrCM5c6cdq2sIMc1HRkUyYMpKrfn0utz91jaVlx8dFkZLS+UfQOvmJK1aZt7KAy/4yi/yde+wORQVKW5iaMNUhB0oqLC1v63frqSwIbYvWCRz4ubfFpt37uOLJV/jwe30SYEckxli+OY0mTNXgpT+9Z3cIqpOrqKrhjhff5/F3P6fO47E7HOUvHfQDaMJUPtVVNaxZvNGWcrvFQPdYISfexfD+nXD6iQP/UrbbP+cu5qYZb7L/QKXdoSg/OXVaiZU0YSrAey9z2i1n2FL2jo272V5QRGH+Lmoqq22JQVnv63VbuPQvM1lbuNvuUJQ/tIWpCVMdcvWdU7nmrqlE2jjdY39RCcP7pzO0dypDB2baFkdQheNETD8VFpdy1dOvMGfpWrtDUapVmjBVAxHh0tum8Ke3buf4c0fbEsOuLXtZ8cUaVi1ch6uzjC7tJNUIlcrqWn79n9n86e351NbpfU2n0i5ZnYepmjFwdG/u/vuPKS4qoXRvOXt27OfeK5/FE8DFbMhx/SgpriAlLZ5t+Tsp2Vvu92s9tXUM6JfBzqIySss67tqk4tYWpj9enL+EdTv28OerziY5PsbucFRTDkxgVtOEqY4qLSuZtKxkuvbODGgu4ZDj+rNq8SYACjfuJiomkmETcsEYtqzbSWlxy8lz1Vfe6Qf9RvaitKzNxTuGSxOm375et4XLnpjJU9dPpX92QE9gUqHg0Baf1TRhqlZtLyjy64kjPXKziU+OxR0ZAcawfvm2w45XH6xhxXfekbgut5DdJ4u66lp2Fxa3+L7h/rSTcLN1bwlXPPkKD18xhcnD+tkdjqqnH0NNmKp1i+auYuiEXLau3U5SagJJ6Ql4jKG6qpaYuCiKd+xn5+Y9xCfHseb7LX69p6fOsHNrMRhD3xE9iYmJpKaqhsoDVWzLb/IM2A6eMI3elmuziqoa/vDax5owHUKXxvPShKla9dY/5lO8qxSA0tIq2HzkY5uGjs9l0+pt9BvSjQ2rtvv/5iIUND7fGPoM70FUVARrbZgXGhIdPOHbxaWji51F/x9rwlQtKy4qbUiWLVm5aCM9+mWxt6idNxtF2Lh6BxjDoHF9qauuIyEtnmFdUvEYQ37+HqqrO9aC7XqZCUx0pF6enERbmDqtRLUiNTOR06Yd59e5WzcUsX9PkEbniLBm6RbWryxk6RfrkLIKVr7xJSkHy8iWavplRNM1JzU4ZYWaXmgCcrCmlsLiktZPVKFnx6IFDvzc6J9wqkUiwi//ciVl+yv4+qPlNgbi/WfXukPdt33HDyI1LR5jDB4DxmPwGAPGUFvroepgDRGRbtwRbjK7pQAG8XXz1f9rjHdf/cAicYl3k/qtUQiN9zX+l0Of7Yb39+0TgbrUSPBYu7B9Z7C3rILzHvk3M2+7nAHddMSs3UTvxWvCVK07UFbJsq/W2xtEM39tFny95qindx3UnczBvTDGUBEbzfr1RSEM7uh69MvgO0/rXdqqeVW1dby+MI+7fzjZ7lCUA1t8VtOEqVq1dulmKssP2h1Gm+xYU8iONYUADD57rG1x6AjZ9nv1q2V0SUng+snjGlrwynp6D1PvYSo/jJw4gLGTh9gaQ3s+q2vmLCYzMzFosShrGQNPvv8lVz39KsXl2rVtC4P3F2H15jCaMFWr3G4XPft3sTUGaUfKNB5Dl6ggBtO20u0quNNZtmkH7y5abXcYYUvXktWEqfywZf1O5sxcaHMU7euKWztnMenp8UGKpQ0c+KHvyN78ZoWu/GQXHSWrCVO1bsF73zvgHmb7Pj111bV0i3UHKZY20FtuQbVhVzHLt+y0O4ywU7/SjxNbmCKySUSWi8j3IrKomeMiIk+JSL6I5IlIwI9i0oSpWuWO6Bz/TdbMXkRqmrWtzIoO/JQVp/r1f2ZTflB/rpay4/5l23oSTjbGjDLGNDfCbwqQ69tuBJ4L9MfQOa6EKqQumn4KZ199vL2JMwijI+uqa8lJsHZgeOl+HaQSbIXFpfx3oY1zglVHcx7wovH6GkgRka6BvJEmTNUqd4Sbmx++hHc3PcHrq/7IpAvsm6bRXmtnLyYlJc6y8jx6vy0k4qJtG8UVtpzaJYv3fs1HIrJYRG5s5nh3YGuj77f59rWZJkzlNxEhITmOW/54KbHx0dYWHqTEU3uwmh7J1l1sPXVGP2Qh8Eneeh59e75OM7GSPYN+MkRkUaOtuYQ40RgzGm/X689E5MQmx5vrngrogqILF6g2i42P5u4Z1zPryQ/YV1RG36Hd6TO4O8npCXTrnUFMfDTJaQneJes8HlZ8swHjMezcupcv3l3Kzi1HPu2kNRX7gvcE6XWzF5H8g+GUlFQG7T1bkhwZzb4avecWTAvXbWHhui3MW7GBe6edynG5Pe0OqdOzaZrHnqPcl2xgjNnu+7dIRN4EjgU+b3TKNqBHo+9zgDY8UukQTZgqIGMmDWbMpMF+ndujf3bD19feNZV1y7Yw56UvWfVdATFxUaRnp7Bo3mpqq2tJzUripKmjKSku54t3l1JbU8ewMb3J+2hJ0GKvqaxmYGoMeZYlzChNmCGydW8JNzz3P66bPJYLjh1G76wOsiB/R2MAj/NuL4hIPOAyxpT5vj4duK/Jae8AN4vIK8BxQIkxZkcg5WnCVJYSEQaO6sXAUb0O279/Txnr87YwfEIuMbHeLtOf/OGHfPPxCmY98N+gx7FuzmISjx1CWVnop8skuCJDXka4+8dni3jp86W88aur6JWpSTMknJcvAboAb/qWTIwAZhpjPhCRmwCMMc8Ds4GzgHygArg20MI0YSpHSMlIZNzkoYftS05L4PRLxtNvUDYv3P0Ky+avoramLijlVR84yICsBPIsSJhRoncxrVBdW8c9r37Mv2+eZnconZITV94xxhQAI5vZ/3yjrw3ws2CUp59k5Xj9Rvbmoffu5KX8p5l4/jiiYtrfYouKjwGPNSujb60qt6SccJeZFM+lE4+4dqpgcfY8TEtowlQdRmqXZO555TZue/aGdr9Xj3H9ydtYHISoWpbVNZnd1XavktT55XbNYOatlzHlmIF2h9JpOXhaiWW0S1a1S12dh5rq2ob7jlY46eLx/Ove1yja2rbRtt1H9iYmOR5cQmxWKpS3fbRuW6X1TAazJ+TlhKuspHhOHNKXa08eS3aqPpEmZA5N8whrrSZMEYnBO0Q32nf+f40x9zY6/jRwrTEmIWRRKseoq/Pw7N2vs2zBWkb8YADffbaS4l0lDDm2H3/490+IS4hp9T02rt5OeUkFmd1SyOyehtvdto6OiMgIpj9+NfdNe8KvhbgTs5LJOeUYVq4vggrfJ9+CZAmwL96A9sgGnQhccOwwrps8Vgf5WMC7lqxmTH9amFXAZGNMuYhEAgtEZI4x5msRGQukhDZE5SS7tuxl9n8WAFC4cXfD/hVf5zN98kNccftZnH7J+MNec7CiivwV20hJT2TTmu08eOMLDcdi46PpP6IH1QdrGX3iQKZcMZEDZZW8/Y/57Nm+j3GnDOXsq47HHXH4wuk/OHcsv3juBh6/6W9HjdUAw84bT8H+Km+ytFhGViJ55XvQFdiD684LJnHK8P5kp2iL0lL6MPTWE6ZvhFH938iRvs2IiBv4E3A5cEHIIlSO0rV3Blf/+hze+ttcSvcdOOxYUeE+nvr1LHZt3UtUdCSL5q1m46pCjDFU+EajpmQcfpGrPFDF8oX5AKxduom3/j6Pmpo6aqtrAVg0dzULP8gjJdaNO9LN/qISeg/N4ep7L+aMH02idN8B3npmDnsK9x0R67ALJ7B8S0kofgx+SRucjqm0PlE7wchPdlJZ6p3nWltdS11to9HNTdcFNgZXhAtPnUHEO/WovufANJn7Fx0bxWV/vhWXS4dfWE1bmH7ew/Qlx8VAf+CvxphvRORW4B1jzA4JwsLYqmMQES679Qz2bN/H7Je+POJ4Xa2HmX/54Kiv37+n5RV7Kg8cOcG/prKKz95e2vD9dx8uo/fQHpx6xQlc/IuzufCWKRQs28yaRRvYtGIrG5dvYfU36yHa4uX7Ghk4tgdrPYf3xR5DLLFLdvkGABrwGOL6plOdGAUCRqAgupadnmrvAEERMAa37/MlBsYtLT08iRjAl2QQvIlEOPJ+k3jfz1NTh6fO21RwuXyP5TaAMd64PJ5DgxN9+3xheN9GGg4d+l6kYb83DmF7YTGetkx093O6UL8R2U4cPNn56T1MwM+EaYypA0aJSAreSaInAhcDk1p7rW/tvxsBevbU5as6ix/fe4H3wrixiKJt+9i+aTcRkW48dZ62XShb4XLJEQ9+7tq3C0PGD2j43u12kTu6D7mj+zTs27tjHwUFuzm1+ACfz13D4m8LAMjumszAwd2Ijokkp0caERFuDlZW859/fhG0C3HupN7Mryjy3szwGeCKIWltCfnLth5+8srCw+sLdGvhvbcEJ8QOqXu/Ltz1wo1tvuetgsGZ0zys1qZRssaY/SIyDzgZb2sz39e6jBORfGNM/2ZeMwOYATB27Fj9iXcSMXHR/PyRSwBva6lkbznJ6Ql8/dFy7rvu6PcV66VnJ7N355HdpSJC7sgeJKXGk9OvCydOHc3gMX248q7zKFy/k/RuqfQf1ZvWejXSu6aS3tU7GOSMKSP59OMV1NV5OOW0YURFH/nfPjE5lr/+5SN/qn6YjIwEuiVHH/rr2xjcO0pJSXGx3xy66ZO1bC/5iza1+f2V18kXH8dPH72cxBRrn2eqDnHiNA+r+TNKNhOo8SXLWOBU4I/GmOxG55Q3lyxVeBCRhnuT408fzpW3T+Glx+Y0e25sfDQDRvbkN8/+iLyv1jP/nSX0G5ZDcmo8dXUexkwaQk6/rCNe12NAN3oMaKntdXTuCBenTxnR4jnnXTiWxd8U8PVX+W167y6psSx/6bMj9mdnJdH/gpEsiqpm9LwiNm/WqSXtUbB8K1HRusSgrbSF6VcLsyvwb999TBfwmjHmvdCGpToqEeGKX57FiB8M4F+PvMuq77xdocPG9+eEc0Zx9tUnNHSpnXTeGE46b4yd4TYQEe685zymX/cCO7bvb/X8iAgXg7snsXLm3GaPlxSVUvL/vuDMs0azvNFoYhWYA6WVlBaXk9k9ze5QVBjzZ5RsHnBMK+foHEx1mOHj+/PYW79gy/qdRES66dY70+6QWuWOcNM9J82vhDmkZwrLZ81v9bxls5cwfPJwVn67MRghhq2+w3posrSTAdFpJbo0ngqtnrnZHSJZAsz89wIW+QYHtaRnj1SWv/J5q+fVW/7ZcoaN79ee0MLelrXbqavTK7atdC1ZTZhK1YuN8295v2RT1+YPc94ny+gzIIuMbroqTSB2bt7D1nUBPcJQBYuxYXMYTZhK+Vxw8bGcec5IXK6jj8BNTIxhxf+OnH/qjw1LCsjumUb3ZgY1qdblLwvnSTX2E2Ms35xGF19XyicmJpLb7zyHa388iVdfXsjHHyzH5RKmnDOKU04fhjvCRUpKLL86uYiCvM0BlZH3yTIABk8cxLaNezngWw1HtWzMKUM5fupou8MIbw5MYFbThKlUE2npCUy/5TSm33Jas8fvf/vXPPqjv7Js/qqAyyhct52oxARNmH5K75pKTJx9KzeFPYOuJYt2ySrVZpk56Tzy4e+4+clrSctu+7MHXG4XmX26sq+oNATRdU5deqbbHUJYE6zvjnVil6wmTKUC4Ha7mPrTM/j78sc4/eqT2vTaYZOGs3FVYesnqgZuXWzdfjpKVrtkVefx4qPvsfDDPA6UHmTqdSdy0fRTj3qux2N4beZCVq3YRlx8NC6Xi5MmD+aYMb2JivL/Y5GQEs8dL0znwlvP4qUH/8eCN75t8fzYhFhWL9Y5mW0xbEIuF958ut1hKAcmMKtpwlSdgjGGbz9dyaY13qkH/3z4XS748clHPEez3oq8rbzw/OGr9Hw8J4/cAdlMu2ICo8f1ISkp1u/y+47oxT2v/pINyzZz74V/omhL80vhRcVHU0WU93HsQE6fDOKi3GxcuZWaqlq/ywsXoycP5f7XbtUF1+2m9zABTZiqkxARMrulsmHFNgD6DO6Gq4WLbN9+WWR3TWbnjsMXgF+/bicP3vsmUVERXHvjJC669Lg2xdFvZC8eeOc33HbiPVQ0GtAjIgw9cyxFO0oob7TofOGmvQBEJyYwYFw229YUeh83LUJVZTXG46G6qrbhkVwdlcslhz3FJiIqouHZly6X4I5wE+f7A8UY0/BMzJv/fIUmS4dw4j1Fq4mx8IcwduxYs2jRIsvKU+GlbN8B5r+9mLikWEafOOiIh1U3taNwH3fePovt2458+HS9p2f8iEFDurc5lnWLC5jzwmd88M+51NXW0WNoT7bvPfJZn/5wR3gThjEGl9uFy/eklvqPrvHN8K7/LIsIgjQcM8aA8XZDYwwez6GEFB0bRUSEC7cYag7WQqMpqMZ3bmPe19aXcUjVwZqGshu/vun3Tc9pTv/hPZh2y+nExEYz7pSh/vyI2k1EFhtjxlpSWAeUHNfNTMi93vJyP8x7wFG/F02YKqzV1Xn4aE4ejz/yfrPH4+KiuO+RaYwc3Sug91+3uICHrnyKjP7dWfn9tvaEGlIRLkNttX8PcQ616Ngo/rPkQRJT4iwrUxNmy5LjupoJ/W1ImMsfbPH3IiI9gBeBbLydxjOMMU82OWcS8DZQP3jgDWPMfYHEo30dynGMMXw/byX/feJ9/nnPa6z5bkPIynK7XUw8YcBRB/pUVFTz73/4v25sUwPG9OWJz+9j8HG5LXYR2+3w9qK9IqOav++sbGRw6ijZWuB2Y8xgYDzwMxEZ0sx5XxhjRvm2gJIlaMJUDvTJywv4zZkP8bc7Z/LKo29z6wn38NPjfsvyL1aHpLyk5Djuuvc8Eo8yyCc2tn3PYUzJTOL6353Po6//nKRUZz4AObKZh2rbZfJFx9odgmqOx4atFcaYHcaYJb6vy4DVQNvvofjJOZ8SpXwOlBw4Yt+GZZu547QHuPWZ6znrhslBL/P4kwYxfmIuWzfvpbKymlUrClm5fCs5PdO5+LLxQSlj6Li+PDTrp9x+/hMN9/ycIrINU2lC7Z0X5lNafIDfPHut3aGoRpw+6EdEeuN9FOU3zRyeICLLgO3AHcaYlYGU4ZxPiVI+p1x+Am/99SN2FOw64tj670M3hzEiwk0f38LoQ4bltHmErD/6Dc3hoVk/4/dXP09F2cGgv3+gjjb9xi4le8vtDkE5Q4aINB74MsMYM6PpSSKSAPwPuM0Y03QJrSVAL2NMuYicBbwF5AYSjCZM5TiJqfE89+1DfPTifL6ZvZTtBUVERrk5dsoxXPKrqXaH125Dxvbhnr/fwF2X/hUrB921pKUntNghItJZCVxh18IFe1objCUikXiT5cvGmDeaHm+cQI0xs0XkWRHJMMY0P1m6BZowlSPFJsRw3k/P4LyfnmF3KCEx8ge5nDbtWD56tbneIxu0MtXDahPOHGl3CKoxA3ic8cddY+Kdo/QCsNoY8/hRzskGdhljjIgci3fszt5AytOEqcJWdVUt5eUHSU6Oa5jraKVr7zyXBbOMap2XAAAIkklEQVSXOapr1gkioyM46fwxdoehDuPMtV2BicBVwHIR+d6377dATwBjzPPARcB0EakFKoFLTYBdO5owVVhasWwrd//qFSoqqklMiuWW289k0inNjUYPnZSMRM6+aiKvP/uppeU2NeiYnuQvC+z5nqFQU1XL879/nZ//8TJHDUYKew5MmMaYBdDynChjzDPAM8EoT/83qrD0zdf5VFRUA1BWWsmjD75DTs80+udmWxrHhDNG2J4wIyLcjlm0oN68Nxex8psNTDx7FBtXFhKbEM1ND1xMWlay3aGFLwcmTKvpPEwVli69YgInn3po2bWa6jr+eP87VFQEtnxdoAaO6klqZstL+IWjmqpatm/czevPfMyiuav44t2lfDhzIdVVzpqOEzbq72FavTmMJkwVluITYrhw2jhc7kO9OZsKdnPHzS+xZXObB88FzOVycaxF66V2dK8+9SHlJRV2hxGmDBiP9ZvDaMJUYWvQkO5MPnXYYfvWr9vJbdNf5PO5oVlVqDnHnDDQsrI6sprqWvY2ebqMspAzl8azlCZMFdauueFEUposV1dWWsn9v3+Dpx//gNra0N/bW71kU8jL6Aw8dR5qqvWZobbQLllAE6YKc9ldU7j5F83P9XznjcW8+vLCkMdQsqcs5GV0BmMmDabv0JAtE6paoy1MTZhKnTR5MA89dildsg8fgelyCT17Z4S8/CFj+4a8jI7uhHOP4f6ZPyMmLtruUMKXJkydVqIUwLjj+vHEc1fz0r8WADD5tGFkZiXRtVtKyMu28rmPHUVsfDSVB7wjlnP6dWH8GSNaffC0CiVnJjCracJUyicjM4nbfnWW5eUmpydYXmZjdqQhd4SL+KQ4ElPiSEiJIyk1noTkOBJT4+g7JIeTfziOPdv3IS4Xmd1SdW1ZuxnA47xRq1bThKmUzXoNyEZEbFuIvT2lRkZHkJQaT2JKPAkpcSQmx5GYFk9SWoL369R4ElJivceTvQkyKS2emLjoVluM3fpktSMyFXTawtSEqZTd0rokc9H0ybat+FOftmLiokhKSyA5LYHE1HiS0rxbYko8ianxJKclePelehNifeJTYUITpiZMpZzgurumcvol4/n0v9/xytMftfv9XC4hKc2b5JLTvVtiSnzD18lpCSSlertC6/dFRUcGoSaqc3LmNA+racJUyiFy+mYxeEzvZo+5XNLQ6kvLTCIlI5GktHhS0hNJTvcmvZR0376MRBKSY3G5dBC8ChIDxoEr71hNE6ZSDpI7ogc/f3gaCclxviTo6wZNS8Dt1gSolJ00YSrlIKmZSZx15US7w1DqSNolqwlTKaWUH3TQjyZMpZRSrTBG52GiCVMppZQ/tIWpCVMppVTrjLYwNWEqpZRqja4lC5owlVJKtab+eZhhThOmUkqp1unCBZowlVJKtcwARluYmjCVUkq1whhtYaIJUymllB+0hakJUymllD+0hYlY+dBaEdkNbA5xMRnAnhCXYSWtj7NpfZzN3/r0MsZkhjqYjkpEPsD7s7TaHmPMmTaU2yxLE6YVRGSRMWas3XEEi9bH2bQ+ztbZ6qPspc8LUkoppfygCVMppZTyQ2dMmDPsDiDItD7OpvVxts5WH2WjTncPUymllAqFztjCVEoppYKuwyZMEblYRFaKiEdExjbaf5qILBaR5b5/Jzfz2ndEZIW1EbesrfURkTgReV9E1vhe94h90R8pkN+PiIzx7c8XkadEROyJvnkt1CldROaKSLmIPNPkNZf56pQnIh+IiB1D85sVYH2iRGSGiKzz/d/7ofWRNy+Q+jQ6x3HXBOU8HTZhAiuAC4HPm+zfA5xrjBkOXAP8p/FBEbkQKLckwrYJpD5/NsYMAo4BJorIFEsi9U8g9XkOuBHI9W2OmX/lc7Q6HQR+D9zReKeIRABPAicbY0YAecDNFsTprzbVx+duoMgYMwAYAswPaYRtE0h9nHxNUA7TYVf6McasBmjaCDHGLG307UogRkSijTFVIpIA/BLvRfk1q2L1RwD1qQDm+s6pFpElQI5F4baqrfUB0oAkY8xC3+teBM4H5lgSsB9aqNMBYIGI9G/yEvFt8SKyF0gC8i0I1S8B1AfgOmCQ7zwPDlrkIJD6OPmaoJynI7cw/fFDYKkxpsr3/f3AY0CFfSG1S9P6ACAiKcC5wKe2RBW4xvXpDmxrdGybb1+HZYypAaYDy4HteFtkL9gaVDv4/p8B3C8iS0TkdRHpYmtQ7dfRrwnKQo5uYYrIJ0B2M4fuNsa83cprhwJ/BE73fT8K6G+M+YWI9A5yqH4JZn0a7Y8AZgFPGWMKghWrP4Jcn+buV1o+hLs9dWrmvSLxJsxjgALgaeAu4IH2xtmGGIJWH7zXixzgS2PML0Xkl8CfgavaGabfgvz7sf2aoDoWRydMY8ypgbxORHKAN4GrjTEbfLsnAGNEZBPeemeJyDxjzKRgxOqPINen3gxgvTHmifbG11ZBrs82Du9SzsHbKrNUoHU6ilG+99wAICKvAXcG8f1bFeT67MXbEnvT9/3rwPVBfP9WBbk+tl8TVMfS6bpkfd1G7wN3GWO+rN9vjHnOGNPNGNMbOB5Y1xE+GEerj+/YA0AycJsdsQWihd/PDqBMRMb7RsdeDbS1BeQ0hcAQEalf1Ps0YLWN8bSL8U7afheY5Nt1CrDKtoDaqaNeE5SNjDEdcgMuwNsqqQJ2AR/69v8OOAB832jLavLa3sAKu+vQnvrgbYEZvBfg+v032F2P9vx+gLF4RzpuAJ7Bt7CGU7aj1cl3bBNQjHe05TZgiG//Tb7fUR7eZJNudz3aWZ9eeEeh5uG9Z97T7nq0pz6NjjvumqCb8zZd6UcppZTyQ6frklVKKaVCQROmUkop5QdNmEoppZQfNGEqpZRSftCEqZRSSvlBE6ZSSinlB02YSimllB80YSqllFJ++P/QAMhEA3Nc6QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1440x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "fig, ax = plt.subplots(figsize=(20, 5))\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"10%\", pad=1)\n",
    "gdf.plot(column='tweetspercapita', ax=ax, legend=True,cax=cax);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}