05.1 - Reproduce Hauer 2016 Inundated Population Results.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import sys, os, time, math, csv\n",
    "import itertools\n",
    "import collections\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "sys.path.append(os.path.join(os.getcwd(),\"simple-maps\"))\n",
    "from simplemaps.SimpleFigures import simpleMap, differenceMap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load counties for which we have migration data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "migration_counties = []\n",
    "f = open(\"data/processed/largestCountyIntersection_2004_2014.txt\")\n",
    "for line in f.read().strip().split(\"\\n\"):\n",
    "    migration_counties.append(line)\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load flooding per block group"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_bg_flood_data(slr_amount):\n",
    "    f = open(\"data/processed/slr_%dft_bg_intersection.csv\" % (slr_amount), \"r\")\n",
    "    header_line = f.readline().strip().split(\",\")\n",
    "    data = {}\n",
    "    for line in f:\n",
    "        line = line.strip()\n",
    "        if line != \"\":\n",
    "            parts = line.split(\",\")\n",
    "            geoid = parts[0]\n",
    "            total_area = int(parts[1])\n",
    "            flooded_area = int(parts[2])\n",
    "            percent_flooded = float(parts[3])\n",
    "            data[geoid] = percent_flooded\n",
    "    f.close()\n",
    "    \n",
    "    return data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "flood_data = [\n",
    "    load_bg_flood_data(i) for i in range(7)\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(flood_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load population per block group"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_bg_population_data(year):\n",
    "    f = open(\"data/processed/regression_estimated_population.csv\", \"r\")\n",
    "    header_line = f.readline().strip().split(\",\")\n",
    "    data = {}\n",
    "    \n",
    "    year_index = None\n",
    "    for i in range(2, len(header_line)):\n",
    "        if int(header_line[i]) == year:\n",
    "            year_index = i\n",
    "            break\n",
    "    if year_index is None:\n",
    "        raise ValueError(\"year must be between 2010 and 2100\")\n",
    "    \n",
    "    for line in f:\n",
    "        line = line.strip()\n",
    "        if line != \"\":\n",
    "            parts = line.split(\",\")\n",
    "            geoid = parts[1]\n",
    "            population = float(parts[year_index])\n",
    "            data[geoid] = population\n",
    "    f.close()\n",
    "\n",
    "    return data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "population_data = load_bg_population_data(2050)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "216330"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(population_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load block group --> county mapping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = open(\"data/processed/regression_estimated_population.csv\", \"r\")\n",
    "header_line = f.readline().strip().split(\",\")\n",
    "bg_to_county = {}\n",
    "for line in f:\n",
    "    line = line.strip()\n",
    "    if line != \"\":\n",
    "        parts = line.split(\",\")\n",
    "        county_id = parts[0]\n",
    "        bg_id = parts[1]\n",
    "        bg_to_county[bg_id] = county_id\n",
    "f.close()\n",
    "\n",
    "block_groups = sorted(bg_to_county.keys())\n",
    "counties = sorted(list(set(bg_to_county.values())))\n",
    "\n",
    "\n",
    "counties_to_index = {county: i for i, county in enumerate(migration_counties)}\n",
    "\n",
    "bg_to_county_index = {}\n",
    "for k,v in bg_to_county.items():\n",
    "    if v in counties_to_index:\n",
    "        bg_to_county_index[k] = counties_to_index[v]\n",
    "    else:\n",
    "        bg_to_county_index[k] = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hauer Results Implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "medium_years = [2055, 2080, 2100]\n",
    "high_years = [2042, 2059, 2071, 2082, 2091, 2100]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "1\n",
      "2\n"
     ]
    }
   ],
   "source": [
    "# medium\n",
    "medium_results_county = []\n",
    "previously_affected_population = {bg_id: 0.0 for bg_id in block_groups}\n",
    "for i in range(len(medium_years)):\n",
    "    print(i)\n",
    "    population_data = load_bg_population_data(medium_years[i])\n",
    "    \n",
    "    affected_population_county = np.zeros((len(migration_counties)))\n",
    "    population_county = np.zeros((len(migration_counties)))\n",
    "    \n",
    "    affected_population = {}\n",
    "    for bg_id in block_groups:\n",
    "        affected_area = 0.0\n",
    "        if flood_data[i][bg_id] != 1.0:\n",
    "            affected_area = (float(flood_data[i+1][bg_id]) - float(flood_data[i][bg_id]) ) / (1.0 - flood_data[i][bg_id])\n",
    "        affected_population[bg_id] = previously_affected_population[bg_id] + ((population_data[bg_id] - previously_affected_population[bg_id]) * affected_area)\n",
    "    \n",
    "    for bg_id in block_groups:\n",
    "        county_index = bg_to_county_index[bg_id]\n",
    "        if county_index is not None:\n",
    "            affected_population_county[county_index] += affected_population[bg_id]\n",
    "            population_county[county_index] += population_data[bg_id]\n",
    "    \n",
    "    medium_results_county.append((population_county,affected_population_county))\n",
    "    previously_affected_population = affected_population"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n"
     ]
    }
   ],
   "source": [
    "# high\n",
    "high_results_county = []\n",
    "previously_affected_population = {bg_id: 0.0 for bg_id in block_groups}\n",
    "for i in range(len(high_years)):\n",
    "    print(i)\n",
    "    population_data = load_bg_population_data(high_years[i])\n",
    "    \n",
    "    affected_population_county = np.zeros((len(migration_counties)))\n",
    "    population_county = np.zeros((len(migration_counties)))\n",
    "    \n",
    "    affected_population = {}\n",
    "    for bg_id in block_groups:\n",
    "        affected_area = 0.0\n",
    "        if flood_data[i][bg_id] != 1.0:\n",
    "            affected_area = (float(flood_data[i+1][bg_id]) - float(flood_data[i][bg_id]) ) / (1.0 - flood_data[i][bg_id])\n",
    "        affected_population[bg_id] = previously_affected_population[bg_id] + ((population_data[bg_id] - previously_affected_population[bg_id]) * affected_area)\n",
    "    \n",
    "    for bg_id in block_groups:\n",
    "        county_index = bg_to_county_index[bg_id]\n",
    "        if county_index is not None:\n",
    "            affected_population_county[county_index] += affected_population[bg_id]\n",
    "            population_county[county_index] += population_data[bg_id]\n",
    "    \n",
    "    high_results_county.append((population_county,affected_population_county))\n",
    "    previously_affected_population = affected_population"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum_medium_results = []\n",
    "sum_high_results = []\n",
    "\n",
    "for i in range(len(medium_years)):\n",
    "    sum_medium_results.append(medium_results_county[i][1].sum())\n",
    "    \n",
    "for i in range(len(high_years)):\n",
    "    sum_high_results.append(high_results_county[i][1].sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEWCAYAAACjYXoKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3XlYlOX6wPHvww6yKijghhtoKceUzB2tXFLTTLNFs2wxy7LFY3XM1Laj7Wm22amsPG36y04u2WJuCZZLimvumoqKooIi28zz++MdkG1ggBlmgPtzXXPJvO8z73sPyNw8u9JaI4QQQlSGm7MDEEIIUf1JMhFCCFFpkkyEEEJUmiQTIYQQlSbJRAghRKVJMhFCCFFpNSaZKKU+VkqdUkptt6Hsm0qpLZbHHqXUuaqIUQghaipVU+aZKKV6AheAz7TWbcvxukeAq7TW9zgsOCGEqOFqTM1Ea70GSC14TCnVQim1XCm1SSm1VinVuoSX3g58WSVBCiFEDeXh7AAcbC4wTmu9Vyl1DfAucG3eSaVUU6AZ8KuT4hNCiBqhxiYTpZQ/0BVYoJTKO+xdpNhtwEKttakqYxNCiJqmxiYTjCa8c1rr9qWUuQ0YX0XxCCFEjVVj+kyK0lqnAQeVUrcAKMM/8s5b+k9CgEQnhSiEEDVGjUkmSqkvMRJDjFLqqFLqXmAkcK9SaiuwAxhS4CW3AV/pmjKcTQghnKjGDA0WQgjhPDWmZiKEEMJ5akQHfGhoqI6KinJ2GEIIUa1s2rTptNY6zB7XqhHJJCoqio0bNzo7DCGEqFaUUoftdS1p5hJCCFFpkkyEEEJUmiQTIYQQlSbJRAghRKVJMhFCCFFpkkyEEEJUWo0YGmyL8+fPc/r0abKzs50dihC1mru7OwEBAdStWxdv76ILeQuHSz0ACXMg6Rs6Rrh1tNdla0UyyczM5OTJkzRq1AhfX18KLEkvhKhCWmtycnJIS0vjyJEjNGnSRBJKVdr7M3wzGkw5YM6x66VrRTNXSkoKYWFh+Pn5SSIRwomUUnh5eREaGkpISAipqallv0jYR+oBI5HkZNg9kUAtSSaZmZn4+/s7OwwhRAGBgYGkp6c7O4zaI2GOUSNxkFqRTHJzc/HwqBUtekJUG56enphMsslplUn6xiE1kjy1IpkA0rwlhIuR38kqln3BoZev0mSilPpYKXVKKbXdyvmRSqkkpdQ2pVRCwZ0RhRBCVIKXY5v6q7pmMg/oX8r5g0C81rod8AIwtyqCEkKIGi92BLh5OuzyVZpMtNZrAKvDN7TWCVrrs5an64FGVRJYNTNv3jyUUuzbt6/YudzcXJRSTJ8+vVj5Q4cOlfted999N40aVfzHsGXLFoYNG5Y/BDQiIoLevXsze/bsQuWUUkyZMsXqdQ4dOoRSKv/h7u5OeHg4I0eO5O+//65wfNXdqlWrUEqxatUqZ4ciXF3Xh8G9hiSTcroX+MHaSaXUWKXURqXUxpSUlCoMC8i6ACv/Da80h+nBxr8r/20cd0EDBw4kMTGRiIiIKr3vhg0b6Ny5M6dPn+aVV17hxx9/5NVXXyUmJoZFixZV6Jr/+te/SExMZOXKlTzxxBN89913DBkyhJwcx3UsurIOHTqQmJhIhw4dnB2KcHV1m8OgNy1P7N9f5ZJDnJRSvTGSSXdrZbTWc7E0g8XFxVXdRvZZF+A/18PZg5CbaRzLOAPrZsHO7+G+X8DbtYYhh4WFERZml83UyuXtt98mODiYn376qdDEtFGjRmE2myt0zebNm9O5c2cAevbsSU5ODlOmTGHTpk35x2sDk8mE1prAwMBa9b5FJZ3aCSijyeuvZYD9/gB2uZqJUioW+A8wRGt9xtnxFJMwu3AiyZObaRxPmF3y65yopGaujIwMHnzwQerVq4e/vz9Dhw4lISEBpRTz5s0rdo0///yTHj164OfnR6tWrXj//ffLvG9qaiohISElznB2c7PPf728v8iPHDlSZtlZs2bRpk0bfH19CQkJIS4urlgNadGiRXTr1g1/f38CAwPp1KkT33//ff753NxcZsyYQevWrfH29iYyMpKJEyeSmXn5/0Nek9wHH3zA1KlTiYiIIDg4mBtvvJGjR48Wut9XX33FtddeS1hYGP7+/lx11VV8+umnxWJXSvHMM88wc+ZMmjVrhpeXF9u2bSuxmUtrzZtvvklMTAxeXl5ERETw8MMPk5aWZtP3VNRQmedh4ydw5VC4eS786yibks2b7HV5l6qZKKWaAN8Cd2qt9zj8hj88DSe2le81f68Hc27J53IzYe0bcGid7dcLbwc3zCxfDBYmk4nc3Nxix2wxduxYFixYwPTp04mLi2PFihWMHDmyxLJpaWnccccdPPbYY0ydOpVPPvmEBx98kJiYGHr37m31Hp06dWLp0qWMGzeOe+65hw4dOth9vk9egmzRokWp5f773/8yceJEpk6dSo8ePbh06RJJSUmFZmC//fbbTJgwgZtuuolPP/0Uf39/Nm/eXCgJjxo1isWLF/PUU0/RtWtXdu3axbPPPsuhQ4f4v//7v0L3nDFjBl27duXjjz/m1KlTTJw4kVGjRhX64D9w4ADDhw/n6aefxs3NjTVr1nDfffdx6dIlxo0bV+h68+bNo3nz5rz22mvUqVOHyMhIzp8/X+y9PvPMM8yYMYPx48dz4403snPnTp599lm2bt3K6tWr7ZbIRTWzaR5kpUG3CQ65fJUmE6XUl0AvIFQpdRSYBngCaK3fB6YC9YB3LWPQc7XWcVUZY5msJZL881XXdt+6desKve6vv/7iiy++YObMmTz55JMA9OnTh4yMDN5+++1i5dPT03n33XfzE0fPnj358ccf+fLLL0tNJpMmTeLPP//kgw8+4IMPPsDX15fu3btzyy23cO+991boQ81sNpObm0t2djbr16/npZdeYtiwYXTsWPp6dYmJicTGxjJ16tT8YwMGDMj/Oi0tjcmTJzN06FC+/fbb/OP9+vXL/3rt2rV8/fXXfPrpp4wePRqA66+/nrp16zJq1Ci2bNlC+/bt88tHRUXxxRdf5D9PSUlh0qRJHD9+nMjISAAmT55c6L316tWL5ORk3nvvvWLJRGvNTz/9hK+vb/6xXbt2FSqTmprK66+/zl133cWcOXPy30NYWBh33nknS5YsYfDgwaV+r0QNlJsN69+DZvEQeZVDblGlyURrfXsZ5+8D7quicCpWI3iludFHYo1fKIxZWvGYymHRokXFRlqZTKYy29B///13tNbccssthY4PHz68xGTi5+dXKGl4e3sTHR1dZtOSr68vixYtYufOnSxbtiy/4/znn39m4cKFLF++vNwT1x544AEeeOCB/OcdOnRg/vz5Zb7u6quv5t133+WRRx5hyJAhdO3aFT8/v/zzCQkJXLhwgbFjx1q9xvLly/Hy8mL48OGFaoR9+/YFYM2aNYWSScFkBdCuXTvAaJLLSyZ79+5l6tSprFmzhhMnTuT3JZXUNNi/f/9CiaQk69evJzs7m1GjRhU6fttttzFmzBhWr14tyaQ22rYA0pNhyByH3cKlmrmqhavvMzrbi/aZAHj4wNX3Vlkobdu2pWXLloWOFW32KklycjIA9evXL3S8QYMGJZYPCQkpdszb27tQP0FprrjiCq644grAWCft/vvvZ/78+SxdupRBgwbZdI08U6ZMYciQIWRkZLBgwQLmzJnDQw89xMcff1zq60aPHk1mZiYfffQR7777Lp6engwYMIA33niDqKgozpwx/kAobRj0qVOnyM7Opk6dOiWez7tGnrp16xZ6npcg8r5vFy5coE+fPvj5+TFz5kxatGiBl5cX7733Xonvx5bReHnNdkXLenh4UK9ePVlYsTYymyHhbWjQFlpc57DbSDIpr64TjFFbRTvhPXwgpJlx3sXlfdCcOnWKZs2a5R8/efKkw+/t4+PDpEmTmD9/Pjt37ix3MmnatClxcUbLZ8+ePUlPT+eTTz5h3LhxdOrUyerrlFL5tZqzZ8/y008/MXHiRG699VZ+//13QkNDATh27Bht27Yt8Rr16tXDx8eHtWvXlng+r7Zhq8TERA4fPszatWvp3v3ywEVrfxDYUovLS2AnTpzgyiuvLHTNM2fOFEtwohbY9zOk7IKbPwQHLmEjPXHl5e1vDP/t9qjRpKWU8W+3R11yWHBJOnXqhFKKBQsWFDpe9Hll5dWAitq9ezdg21/aZZk5cya+vr4899xzNr8mJCSEW2+9lREjRrB9u7GyT9euXfH392fuXOuLLvTv35/MzEzOnz9PXFxcsUd5k0lGRgZgLHiY5+zZs/zvf/8r13UK6ty5M15eXnz11VeFjn/99dfk5ubSq1evCl9bVFPrZkFgI2MUlwNJzaQivP2h92TjUQ21bt2aO+64g2effRaz2UzHjh359ddfWbx4MWC/Ybtjx44lLS2NYcOG0bZtW0wmExs2bOCVV16hRYsWDB1a+D/37t27WbhwYbHrXHed9ap5eHg448eP57XXXmPTpk1WO+LHjh1LQEAAXbp0oX79+uzZs4fPP/88v78jICCAGTNm8MgjjzBs2DBGjhxJQEAAW7ZswcfHh0ceeYRevXpx++23M3z4cJ544gk6deqEm5sbhw4dYtmyZbz88stER0fb/P3p2rUrgYGBjB8/nueee46LFy/y4osvEhoaWuIoLVvUrVuXiRMnMmPGDOrUqcOAAQPYtWsXU6ZMoXv37gwcOLBC1xXV1NGNcHgd9Jvh0NnvgDFCpLo/OnbsqEuzc+fOUs9XN5988okG9N69e4udy8nJ0YCeNm1asfIHDx7MP3bx4kU9btw4HRISouvUqaNvvPFGvWTJEg3o7777Lr/cXXfdpRs2bFjsPvHx8To+Pr7UOJcvX65Hjx6to6Ojtb+/v/by8tLNmzfX48eP1ydOnChUFrD62LBhgz548KAG9IcffljsPikpKTogIEAPHjzYaizz5s3T8fHxOiwsTHt5eemoqCj92GOP6fPnzxcqt2DBAt2pUyft4+OjAwICdKdOnfTixYvzz5tMJv3WW2/p2NhY7e3trQMDA3VsbKyeNGmSPnfunNZaW4115cqVGtArV67MP7ZixQrdvn177ePjo5s3b65nzZqlp02bpo1fzcLfn2eeeabY+yrpmmazWb/xxhs6Ojpae3p66vDwcP3QQw8Ve6+uoKb9brqcr0ZpPaOx1pnpJZ4GNmo7fQ4r43rVW1xcnN64caPV87t27aJNmzZVGFH19Nprr/Hkk09y6NAhmjRp4uxwRC0gv5sOdGY/vN0RejwB100tsYhSapO20/QLaeaqpZYsWcL27dtp3749bm5urF27ltdee40RI0ZIIhGiJkicA+5e0OmBssvagSSTWiogIIDvvvuOmTNncvHiRRo2bMiECRPK1ZEthHBRF07Bn/+Ff9wGASUP+bc3SSa1VHx8POvXr3d2GEIIR/hjLpiyoesjVXZLGRoshBA1SdYF+ONDaD0QQltV2W0lmQghRE3y53zIPGfMfatCNiUTpZSXUmqaUmq3UipDKWUq8ih7DQ8hhBCOZcqFxHegSRdobH1FCEewtc/kVWA8xs6H3wJZDotICCFExez8Ds4fgRtervJb25pMhgPTtNYvOTIYIYQQFaQ1rHsLQqMhun+V397WPhN/INGRgQghhKiEA6uMzf66TgAnbIBm6x0XAz0dGYgQQohKWDcL/MON/d2dwNZk8jZwu1JqqlIqTinVvOjDkUGKwvL2dFdKsWdP8d2NV69enX/+l19+sdt9lVJMnz49//n06dPLvbmVI508eZIJEyYQHR2Nr68voaGhdOzYkUcffZSsrMvdfL169Sq05HtJoqKi8r+HSimCg4Pp06cPv/32m6Pfhksr+n9AuIjkrXBgJXQeBx7FN1arCrYmk0SgFTAd+B3YW8Kj1jmVlsmIDxI5lW7bJlH2FhAQwOeff17s+KeffkpAQIDD73/fffeRmOgarZ9paWlcc801LFmyhMcff5xly5bxwQcfMGDAABYvXsylS5fKfc1+/fqRmJjIunXrmDVrFvv372fAgAGF9oSvbRITE7nvvqrbDFXYKOFt8AqAjmOcFoKtHfD3YKzgKgqYvWIvGw6lMnvFPl68qeQNlRzp5ptvZv78+Tz//PP5NYRLly6xcOFChg0bxrx58xx6/0aNGpW6M2FVWrhwIYcPH2bLli384x//yD8+bNgwnn/++QpdMzQ0NH8L5K5du9KyZUu6d+/OV199xdNPP22XuKuLrKwsvL29y9wSWjjB2cOw/Vvo/CD4BjstDJtqJlrreVrrT0t7ODpQV3MqLZMFm46iNSzc+LdTaid33nknhw8fLtT0smjRIsxmM8OGDSvxNatXr+a6664jICCAOnXq0K9fv/wNovKYTCamTJlCREQEfn5+9OrVix07dhS7VtFmrkOHDqGUKpbEVq1ahVKKVatW5R/La2pavnw57du3x9fXl6uuuorff/+d3NxcJk+eTEREBHXr1uXuu+/m4sWLpX4v8rajDQ8PL3Yur6mqsjp06AAYe7iX5ccff6Rr164EBQXh7+9PTExMsaS2detWhg4dSr169fD19SUmJoYZM2YUKvPtt9/SuXNn/Pz8CA4O5pZbbil2/6ioKEaNGsVXX31FmzZtqFOnDnFxccWa5DZs2MDw4cNp1KhR/v0mT55crNaW97NZvHgxV111Fd7e3rz77rtAyc1cy5cvp0uXLvj6+hIUFMRNN93EX3/9Veb3SNjJ+neNTfo6P+TUMMq1NpcyfiOvAOoCqcBOXY3XsH9u8Q52Hk+r0GsPnr5Idq4ZgKxcM4Nm/0az0JL3Bi/NFZGBTLvxyrILlqBp06b07NmTzz//nB49egDw2WefMXToUPz9i+/4uHTpUoYMGcLAgQOZP38+AC+//DI9evQgKSmJxo0bA0aS+Pe//80TTzxB37592bhxI4MHD65QjKXZt28fkyZN4plnnsHf358nn3ySwYMHM3jwYHJzc5k3bx67du1i0qRJ1K9fn1deecXqtfK27L3tttt4+umn6d69u9W92isqr3mrRYsWpZY7cOAAgwcPZvjw4UydOhUvLy/27t3LgQMH8sv88ccf9OrVi5YtW/Lmm2/SqFEj9u7dS1JSUn6Z999/nwcffJAxY8YwdepU0tPTmT59OvHx8SQlJRVqyly7di1//fUXL7zwAj4+Pjz77LMMGjSIQ4cOERxs/LV65MgR2rdvz913301AQAA7duzg+eef58CBA8V2ZtyzZw8TJkzg2WefpXnz5la3+12+fDkDBw7k2muv5euvv+bChQtMnTqV7t27s2XLFho2bFiu77Eop4xU2PwZtLsFgpz7vbY5mSil7gNeBMIKHD6llJqitf7I7pG5sOxcMykXsvLb/TSQciGLhsG+eHlU7ZC80aNHM3HiRGbPns3Zs2f55Zdf+OGHH0os++ijjxIfH19oW9jevXvTvHlzXn/9dd566y3Onj3Lm2++ydixY3nttdcA6Nu3L+7u7nZv2jlz5gwJCQk0b26M3zCbzQwZMoSDBw/mDxzo168fa9asYcGCBaUmk549e/Lcc8/x4osv0r9/f9zd3Wnfvj2DBg3isccey/9ALQ+tNbm5uZjNZvbt28eDDz5Iq1atuOeee0p93ebNm8nOzua9994jMDAQgGuvvbZQmX/+85/Uq1eP9evX4+fnV6zMhQsXeOqppxgzZgwff/xx/vFOnToRExPDRx99xGOPPZZ/PC0tjS1bthASEgIYNbSrr76aZcuWcccddwAUqq1qrenWrRuBgYGMHj2ad955h3r16uWfP336ND/99BPt27cv9b1OmTKF5s2b88MPP+DhYXycdOnShejoaF5//XXeeOONUl8vKmnDR5CTUaULOlpjUzJRSo0E5gIrgPnACSAcGAnMVUplaK2/dFiUDlLRGsGURdvYfvw8OabLlTIPN8WVDYOqvO/klltu4eGHH2bx4sUcPnyY8PBwrrvuOtasWVOo3N69e9m/fz+TJ08mN/fy6jd+fn506dIlv/y2bdu4ePEiI0YUHl6Y9xe/PUVHR+cnEjC2EwYjgRTUunVrFi9ebOzmVkpz1dSpUxk7dixLly7lt99+Y9WqVTz33HP85z//YdOmTTRoUL6luL/44gu++OKL/Od169YlISEh/wPbmvbt2+Pp6cltt93GPffcQ8+ePalfv37++YyMDNatW8ekSZPyE0lRiYmJpKWlMXLkyEI/r8aNG9O6dWvWrFlTKJl06dKlUFzt2rUDCjfJpaWl8dJLL7Fw4UL+/vtvcnJy8s/t3bu3UDKJiooqM5FcvHiRzZs3M3ny5PxEAtCsWTO6devG6tWrS329qKScS/DHB9CqLzSo2GeZPdn6Z/STwH+11n0sfSQ/Wv7tC3wBPOW4EF1LXl9JwUQCkGPSTuk7CQgI4KabbuLzzz/ns88+Y+TIkSXu4X7q1CkA7r33Xjw9PQs9lixZwpkzZwBITk4GKPbBW94PYlsU/VD28vKyejw3NxeTyVTmNcPDw7n33nv55JNPOHjwIHPmzOHYsWO8+uqr5Y7vhhtuYMOGDSQkJPDWW29x6dIlbr75ZjIzS/8Zt2zZkh9//BGz2cydd95JeHg4nTt3zv9wPXv2LGazudTBC3k/r+uvv77Yz2vbtm35P688RZuhvL2N4aEFYx0zZgzvv/8+EyZM4Oeff2bDhg288847xcoBRERElPoe896H1rrEsuHh4fn9WMJBtn4JF1OMSYouwNZmrhiMhFKS+cB3tlxEKfUxMAg4pbUu9ie8pU9mFjAAyADu1lpvtjHGKjF7xV7MVrqJTFo7ZWTX6NGjGThwIGazmS+/LLmCmPdX54wZM7j++uuLnc/7IM/7YDh58iRXXnn5r52TJ0+WGYePjw8A2dnZhY4X/eCrSuPHj+fZZ59l586d5X5t3bp1iYszdjTt0qULQUFBjBkzhrfffptJkyaV+trevXvTu3dvsrKyWLduHVOnTmXgwIEcOnSIkJAQ3NzcOHbsmNXX5/285s2bV+jnkKe8Q78zMzP53//+x/Tp03n00curyW7btq3E8rYMWAgJCUEpxYkTJ4qdO3HihNV+FmEHZpMxHDiyA0SVPmeqqthaM0kHrP0Z1chy3hbzgNIWjbkBYz5LK2As8J6N160ym4+cK1YryZNj0mw+fLaKI4I+ffowYsQIxo0bV+IHD0BMTAxRUVHs2LGDuLi4Yo/Y2FgAYmNjqVOnDt98802h1xftoC1JgwYN8Pb2LjY6bOnSpRV8Z7Y7efIkZrO52PHk5GTOnz9v01/aZbnrrrvo0KEDr776KhkZGTa9xtvbm2uvvZYnn3ySixcvcvDgQfz8/OjevTvz58+3Ov+la9euBAQEsG/fvhJ/XjExMeWKPSsrC5PJhKenZ6HjlRk+XqdOHTp27MiCBQsK1RoPHz5MQkICvXr1qvC1RRl2L4XUA8Yy8y4ycdjWmskPwL+VUnu01mvzDiqlumB0ypfc41uE1nqNUiqqlCJDgM8sI8TWK6WClVIRWutkG+N0uGWP9nB2CMW4u7tbrZHkUUrxzjvvMGTIELKzsxkxYgShoaGcPHmShIQEmjRpwhNPPEFwcDCPP/44L730EgEBAfTt25cNGzbw0Udlj7FQSnHrrbfy0UcfER0dTUxMDEuXLi00JNhRPv/8c+bOncvIkSPp1KkTfn5+7Nmzh9dffx0vLy/Gjx9fqPyZM2dYuHBhsevExsYSHR1d4j2UUjz//PMMGjSI9957j4kTJ5ZY7v3332fNmjUMGDCAxo0bc/r0aWbMmEFkZCRt2xq11tdee434+Hi6dOnCxIkTadSoEQcOHGDLli28/fbbBAYG8uqrrzJ+/HhSUlK44YYbCAoK4tixY6xevZpevXrld6zbIigoiM6dO/P6668TERFBaGgoH3/8cam1I1u88MILDBw4kEGDBvHQQw9x4cIFpk2bRlBQkNXvj6ikvAUdQ5pBmxudHc1lWusyHxid7XsAE3AEYxb8Ycvzv4AGtlzHcq0oYLuVc0uA7gWerwDiyrpmx44ddWl27txZ6vnq5pNPPtGA3rt3r9UyK1eu1ID++eefCx1PSEjQAwcO1MHBwdrb21s3bdpU33rrrTohISG/TG5urn7mmWd0gwYNtI+Pj46Pj9c7duzQgJ42bVp+uWnTpmnjv9BlZ8+e1aNGjdL16tXTISEh+oEHHtBLlizRgF65cmV+ufj4eN2tW7dCrz148KAG9IcffljoeN59cnJyrL7fnTt36scee0y3b99e161bV3t4eOjw8HA9bNgwvWnTpkJl4+PjNcYgvGKPV199VWutddOmTfXIkSNLvFeXLl10eHi4zsjIKPF8QkKCHjx4sG7UqJH28vLS4eHhevjw4Xr37t2Fym3evFkPGjRIBwUFaR8fHx0TE6NnzpxZqMzSpUt1r169dEBAgPb19dUtW7bUY8aM0Tt27MgvYy3Woj+vgwcP6v79+2t/f38dFhamx48fb/PPxto1tdb6hx9+0J07d9Y+Pj46MDBQDx48uNh7taam/W5WiYO/aT0tUOs/Piy7bBmAjdrGz+6yHkrbOE1EKeWHMRO+B5fnmawG5mmtbavzG9eJApbokvtMlgAztda/WZ6vAJ7SWm8soexYjKYwmjRp0vHw4cNW77lr1y7atGlja4hCiCoiv5sV8MWtcHQDPL4DPH0rdSml1CatdZw9wrJ5noklYcyxPBzlGNC4wPNGlmMlxTMXY7gycXFx1XbipBBC2OzULtizHHpNrnQisTdX2wP+e2C0MnQGzmsX6i8RQginSngbPP2g0/3OjqQYqzUTpdQBYKjWeqtS6iClL/SotdalrzFhXPNLoBcQqpQ6CkwDPC0XeB9YhjEseB/G0GDnLYEphBCuJO04JH0DcWPAz/WGXZfWzLUaSCvwdaWbkrTWt5dxXmPsNS+EEKKg9e+BNkEX1/yItJpMtNZjCnx9d5VE40C6jKU4hBBVy9bBPwLIPA+b5sGVQyEkytnRlMimPhPLDouRVs5FKKWm2jcs+/L09KzQ5khCCMe5dOlS/rIvogyb5kFWmsssnVISWzvgp2F9Bnyk5bzLql+/PseOHSMjI0P+GhLCibTW5OTkkJq0xPJwAAAgAElEQVSaytGjRwstLimsyM02mriaxUNk6YtvOpOtQ4NLax8KAbJKOe90ecuAHz9+vNBKqUKIqufh4YGPjw9NmjTJX89NlGLbAkhPhiGOnJVReaWN5uoFFNyE4QGl1KAixXyBgUDxbfhcTGBgYH5SEUKIasFshoTZ0KAttLjO2dGUqrSaSTwwxfK1puRhutnATsB1G/KEEKK62vsTpOyGmz90mQUdrbHaZ6K1fk5r7aa1dsNo5uqc97zAw0dr3UFrnVh1IQshRC2RMBsCGxmjuFycTX0mloQihBCiqvy9AQ6vg34zwN2z7PJOZvPaXHmUUvWBYr1mWusjJRQXQghREQmzwCcYOox2diQ2sXUPeDeMfUseAIKtFHO3V1BCCFGrndkPu5ZAjyfA29/Z0djE1uarxzCWOXkdo//k3xjJ5SCwH3C9VceEEKK6Sngb3L2g0wPOjsRmtiaTMcDzwMuW54u01tOANhhLxDdxQGxCCFH7XDgFW76Af9wGAQ2cHY3NbE0mzTF25DIBuRjzS9Ba5wBvYWyaJYQQorL+mAumbOj6iLMjKRdbk8l5Lne6HwdiCpzzwNh5UQghRGVkXYA/PoTWAyG0lbOjKRdbR3P9CVwB/Gh5PKeUuoRRS3kJ2OyY8IQQohb5cz5knoNujzo7knKzNZm8hdHUBcaijh2A/1qeHwYetnNcQghRu5hyIHEONOkCjTs5O5pys3XS4s8Fvj6hlOoEtAD8gF2WvhMhhBAVteM7OP83DHjV2ZFUSLknLUL+joj77ByLEELUTlobkxRDo6FVP2dHUyGlrRrcszwX0lqvqXw4QghRCx1YCSe2weA54FY9V68qrWayCtv2fVeWcjIDXgghKmLdLPAPh9gRzo6kwkpLJr2rLAohhKitkrfCgVVw/XTwqL7bGFtNJlrr1VUZiBBC1ErrZoNXAHQsacuo6qN6Ns4JIURNcPYw7FgEHe8CX2tr6FYPtq4a/GsZRbTW2rX3lBRCCFez/l1jB8XODzk7kkqzdWiwG8U74+thLKuSAuyxZ1BCCFEjpR6AhDmQ9A1kXwA01GsJpixnR1Zptk5a7FXScaVUC+A7jCXpbaKU6g/Mwhj99R+t9cwi55sAn2Lsm+IOPK21Xmbr9YUQwiXt/Rm+GW3MdDcXmOedegje6wYjPoNWfZwWXmVVqs9Ea70fmAnYNGVTKeUOvAPcgLHW1+1KqSuKFJsCfKO1vgq4DXi3MjEKIYTTpR4wEklORuFEAqBzjePfjDbKVVP26IBPAaJtLNsJ2Ke1PqC1zga+AoYUKaOBQMvXQRirFAshRPWVMMeokZTGlAOJ71RNPA5QqWSilKoHPIGx26ItGgJ/F3h+1HKsoOnAKKXUUWAZUOKi/kqpsUqpjUqpjSkpKeWKWwghqlTSN8VrJEWZcyDp66qJxwFsHc11kOId8F5A3jZgw+wY0+3APK3160qpLsDnSqm2WmtzwUJa67nAXIC4uDhbZuoLIYRzZF+wbzkXZOtortUUTyaZGMvPL7D0ndjiGNC4wPNGlmMF3Qv0B9BaJyqlfIBQ4JSN9xBCCNfi5Q/Z6baVq6ZsHc11t53utwFopZRqhpFEbgPuKFLmCHAdME8p1QZjh0dpxxJCVF+xI2Dzp2DOtV7GzRNib626mOysQn0mSqmwirxOa52LsZHWj8AujFFbO5RSzyulBluKTQTuV0ptBb4E7rYseS+EENVTVxv2D3T3hC7jHR+Lg9i8n4lSKh54HmNElpdSKhv4HZhanuXnLXNGlhU5NrXA1zuBbrZeTwghXN6FFDCbQLmDcivcGe/maSSSEZ9B3ebWr+HibKqZKKVuAX4F6mPMKZkAvIbRAf+rUmq4wyIUQojqLOsCLBoLwY1h7CpjHS7vAGMZFe8A4/mD66r1hEUAZUsLklJqF7AXuKngqCqllBvwPdBCa93GYVGWIS4uTm/cuNFZtxdCCOu+fwQ2fw5jlkHTrs6OphCl1CatdZw9rmVrn0kz4L0ShueaMWaoR9kjGCGEqFF2L4PNn0G3R10ukdibrclkL2Ct0z0M2Q9eCCEKu5Bi1EoatIPek50djcPZmkyeAZ5TSl1d8KBS6hqMGev/snNcQghRfWltJJKsdLh5brXeQdFWto7mmoQx32O9Uupv4CRG53tjy9dPKqWetJTVWut4u0cqhBDVxebPYM8P0O/f0KDoWrY1k63JxATstjzyHLQ8hBBC5DmzH5b/C5r1hGsedHY0VaZS+5kIIYQowJQLi8aBuwfc9B641Z6d0W2etCiEEKIMv70JR/+AYR9BUCNnR1OlbE6bSqkIpdRrSqkNSqn9ln9fUUqFOzJAIYSoFo5thtUzoe0waFf75nHbOgM+GtiCMfP9AvCH5d9HgS1KqVYOi1AIIVxddgZ8Oxbq1IeBrzs7GqewtZnrZSANuEZrfSjvoFKqKfCT5fzNdo9OCCGqg1+mwZm9cOd34Bvi7GicwtZmrt7AswUTCYDW+jDGPJPe9g1LCCGqiX2/wB9zjZFbLWrvR6GtycQLsLazS7rlvBBC1C4ZqfDdeAhrDddPc3Y0TmVrMtkCPGJZ2DGfUkoBD1nOCyFE7aE1LHkMMs4Ys9w9fZ0dkVPZ2mfyPLAE2KWU+hpIBsKBW4BWwEDHhCeEEC4q6WvY+T+4bhpE/MPZ0TidrZMWlyulBgEvYqzTpTD2hN8EDNJa/+S4EIUQwsWcOwLLJkGTLsaKwML2SYta6+XAcqWUHxACnNVaZzgsMiGEcEVmEyx6ELQZhr4Pbu7OjsgllJpMlFLewA1AS+AcsERrfQKQJCKEqJ0S34HDv8GQdyAkytnRuAyryUQp1RD4BYjGaNYCuKiUukFr/VtVBCeEEC7lxHb49QVoPQjaj3R2NC6ltNFcL2BsfDUGuBK4ETgKvF0FcQkhhGvJyTRmufsEw42zjD3cRb7Smrl6AVO11p9Znu9SSp0B1imlQrTWZx0enRBCuIqVL8KpHXDHN1An1NnRuJzSaiaNKT5/5E+MJq/GDotICCFczcG1kDAH4u6B6H7OjsYllZZM3IHcIsfynsvS9UKI2iHzvLFHSd3m0PdFZ0fjsspKCmMt80vy5M0veVAplVzguNZa27SWgFKqPzALI1n9R2s9s4QyIzDW/NLAVq31HbZcWwgh7G7ZJEhPhnt/Aq86zo7GZZWVTO6xcvzeIs81UGYyUUq5A+8AfTA68zcopb7XWu8sUKYV8C+gm9b6rFKqflnXFUIIh9j+rTHTPf5paBTn7GhcmtVkorV2xH6TnYB9WusDAEqpr4AhwM4CZe4H3snr4Ndan3JAHEIIUbq047DkcWjYEXr+09nRuLyq3qC4IfB3gedHLccKigailVLrlFLrLc1ixSilxiqlNiqlNqakpDgoXCFErWQ2w3cPgSkbhs4Fd09nR+TyXHG3ew+MxSN7AbcDHyqlgosW0lrP1VrHaa3jwsLCqjhEIUSNtuFDOLDS6HAPbensaKqFqk4mxyg8rLiR5VhBR4HvtdY5WuuDwB6M5CKEEI6X8hf8PBVa9jGGAgubVHUy2QC0Uko1U0p5AbcB3xcp8x1GrQSlVChGs9eBqgxSCFFL5WbDt/eDpx8MmSOz3MuhSpOJ1joXeBj4EdgFfKO13qGUel4pNdhS7EfgjFJqJ7ASmKS1PlOVcQohaqnVL0PyVmO5lIBwZ0dTrSittbNjqLS4uDi9ceNGZ4chhKjOjvwOn/SHf9wON73r7GiqhFJqk9baLmOeXbEDXgghqlZWOiwaC0GNoH+xedTCBqUtQf9rOa6jtdbX2SEeIYSoej9OhrOHYcwy8Al0djTVUmkz4N0wZrbnicHY9/0QcBJoAERh7Af/l2PCE0IIB9u9DDZ/Bt0fh6ZdnR1NtVXaDPheeV8rpW7CWE+ri9b69wLHrwG+tpwTQojq5cIp+P4RCG8HvSY7O5pqzdY+kxeAZwsmEgDL8+mALKUphKhetDYSSVY63PwheHg5O6Jqzdal5FsB1tYsOYWxR7wQQrim1APGfiRJ30D2BfDyh4hYOLwO+s2A+m2cHWG1Z2syOQg8APxQwrkHMPpRhBDC9ez9Gb4ZDaYcMOcYx7LTjUSi3KBuC+fGV0PYmkyeA/6rlNoOLORyB/xwoDUw0jHhCSFEJaQeMBJJTkbJ57UZFt4ND64zNr8SFWZTn4nW+iugH3AeY6+Rdyz/ngP6aa2/dliEQghRUQlzjBpJaUw5kPhO1cRTg9k8aVFr/YvWuhvgizFE2Fdr3V1rvcJh0QkhRGUkfXO5acsac46xAZaolIrMgPfDSCjudo5FCCHsK/uCfcsJq2xOJkqpQUqpzRhNXQeAdpbj/1FKyR7tQgjX4+Vv33LCKpuSiWXS4v+A08BTQMF1mQ8Cd9k/NCGEqKTofmWXcfOE2FsdH0sNZ2vNZBrwida6L/BWkXPbgbZ2jUoIISpr61ewe2nZ5dw9oct4x8dTw9maTNpgLJsChdfrAjgL1LNbREIIURlZF2DRg7DoAYi8Cm5639jsyq3IPu5unsbxEZ/JsGA7sHWeSRoQauVcFNZnxwshRNU5sQ0WjIEz+yD+Kej5JLh7QJNrjOG/SV9fngEfe6tRI5FEYhe2JpOfgX8ppX4A0i3HtFLKG2PnxJJmxgshRNXQGjb8B358BnxD4K7voVnPy+frNoeBrxsP4RC2JpNngD8wlppfhtHU9TQQCwQBNzkkOiGEKMuls/C/h2H3EmjZB4a+D3WsNaQIR7F1BvwhoAOwBOgDmICewHrgGq31cUcFKIQQVh35Hd7vAXuWQ98X4Y5vJJE4ia01E7TWR4F7HRiLEELYxmyGdW/Cry8ZW+3e8xM06ujsqGo1W+eZ/KqUam3lXHQ5t/gVQoiKSz8J84fCiufhisEwbq0kEhdga82kF2BtY+QAIN4u0QghRGn2rTCG/GZdgBtnQ4fRoFTZrxMOZ3MzF8Xnl+RpAcjCNkIIxzHlwMqX4Lc3Iaw13LVYNrRyMVaTiVJqDDDG8lQDc5VS6UWK+WLMfpeVg4UQjnHuCCy8F47+AR3ugv4zwcvP2VGJIkqrmZgxRm2BsRZXwed5zgDvAS/bPzQhRK2383v4/mFjHsnwj6HtMGdHJKywmky01p8CnwIopVYCD2mtd1X2hkqp/sAsjCXs/6O1nmml3DCMXR2v1lpvrOx9hRDVSE4m/PSMMRExsoORSOo2c3ZUohQ29ZlorXvb42ZKKXeMXRr7AEeBDUqp77XWO4uUCwAeBX63x32FENVIyh5YOAZObocuD8N108DDy9lR1Tin0jLxqNc4xl7Xs3Vo8JtKqc+tnPtcKfWqjffrBOzTWh/QWmcDXwFDSij3AkbTWaaN1xVCVHdaw5//hbnxkJ4MdyyAfi9JInGQ2Sv24ubhbbeNXGxdNXgw8JOVcz9i+3IqDYG/Czw/ajmWTynVAWistS517Wil1Fil1Eal1MaUFFlnUohqLSvdGPL7v4egYUcY9xtE93V2VDXWqbRMFmw6WnhnqkqydWhwQ+CIlXPFEkJFKaXcgDeAu8sqq7WeC8wFiIuLszZsWQjh6o5vMZq1zh6C3s9Aj4ngJruCO9LsFXsxa/t+bNpaMzkLtLRyriW2zzM5BjQu8LyR5VieAIyhxquUUoeAzsD3Sqk4G68vhKgutIb178NHfYwO97uWQPyTkkgcLK9WkmOybzKxtWbyCzBFKbVEa30y76BSqgEwGWOJeltsAFoppZphJJHbgPz947XW5ymwb4pSahXwTxnNJUQNk5EK/xsPfy2D6P4w5F2oI3vsOdrfqRlM+PJPsnLNdr+2rcnkWYxEsFcptYTLTVuDMDrJp9hyEa11rlLqYYx+FnfgY631DqXU88BGrfX35X0DQohq5nAC/N99cOGUMQHxmnGyJIoDHTt3iWVJySxJOs7Wo+cddh9bhwYfUkpdDTyPMay3HnAaWARM01oftvWGWutlGHuiFDw21UrZXrZeVwjh4swmWPsGrPo3BDeF+342ttUVdnf83CWWbUtm6bZk/jxyDoDYRkFc1TiY7cfP272JC8q3BP0hYLTdIxBC1HxpyfDt/XBoLbS7BQa+AT7W1o4VFXEyLdNIIEnJbDx8FoArIwN5sn8MA9tF4OvpTo9XVjokkUD5FnrMG211BUbNZKPW+qJDohJC1Bx7fzGG/WZfhMFz4KpR0qxlJ6fSM1m+/QRLtiaz4XAqWkPr8AD+2TeagbGRNAutk192yqJtdh/BVZDNyUQpNR6YhtFBroGrgc1Kqe+AX7XWsx0TohDCpaUegIQ5kPQNZF8AL3+IHQGdHoAt8yFhNtS/Em75BMLsNuG61jp9IYsftp9gadJxfj9oJJCYBgE8fn00A9pF0LJ+yfMQNx8557BaCYDSNmQqpdT9GAs6fowxefEbIE5rvVkpNREYrLV22p4mcXFxeuNGGfAlRJXb+zN8M9pYIt6cc/m4mwdos/GIu9eYye7p67w4q7nUi9ks336CpduOk7j/DGYNLcLqMCg2koGxEUQ3CKjQdZVSm7TWdpl6YWvN5Angda31U5b1tQraDUyyRzBCiGok9YCRSHIyip8z5xr/untB14clkVTA2YvZ/LTzBEuSkknYfwaTWdMstA7je7dkUGwk0Q38US7UXGhrMmmGMZy3JBeBYPuEI4SoNhLmGDWS0mgNie/AwNerJqZq7nxGDj/uPMHSpGTW7TtNrlnTtJ4f4+KbM7BdJG0iAlwqgRRkazI5DURZORdD4VnsQojaIOmbwk1bJTHnQNLXkkxKkZaZw887TrJ0WzJr96aQY9I0ruvLfT2aMyg2gisjA102gRRkazJZAky1zEjPm1OilVKhwOPAdw6ITQjhyrJtXEXJ1nK1SHpmDr/sOsnSpGTW7DlNtslMw2BfxnRrxsB2EcQ2CqoWCaQgW5PJFKA3sB1jjxENzAZaA6cwJjMKIWqDlD2QMAvjY8AGXnZb5bxau5CVywpLAlm1J4XsXDMRQT6M7tKUgbERtG8cXO0SSEG2zoA/bVls8TGgH7Df8to5wJta6zTHhSiEcAl//wG/vQV/LQUPHwhrA2f2Xu5sL4mbJ8TeWnUxupiM7Fx+3X2KpUnJ/Lr7FFm5ZhoEejPymiYMio3gqsYhuLlV3wRSkNVkopR6AyNR/K2U6gls1lq/gLFxlRCiNjCbYe9PsO4tOJIIPsHQ80m45gHISoP3upWeTNw9ocv4qovXBVzKNrHqr1MsSUpmxe6TZOaYCQvw5rarGzMwNpK4pjUngRRUWs3kMYydEP8GVgJdgD+qIighhJPlZsO2BcaEw5TdENTYWJTxqjshb3O+OqEw4jMr80w8jUQy4jOo29w576EKZeaYWPVXCku3JbNi10kysk2E+nsxvGMjBsVGcnVUXdxrYAIpqLRkkgK0wkggCpsbSIUQ1VZWOmyaB4nvQvpxY+b6zR/ClUON5FBUqz7w4Dpj+G/S1wVmwN9q1EhqcCLJyjWxZs9pliYd5+edJ7mYbaJuHS9uuqohg9pF0KlZXTzcbd0yqvqzOgPesuf7rcBxoAmQDGRZuY7WWrdwSIQ2kBnwQlRS+kn4/X3Y8BFknYeoHtDtMWh5nayjVUB2rpnf9qWwZGsyP+88SXpWLsF+nvS/MpyBsRF0aV6vWiWQqpoB/wCwDWiDsVpwEsZ8EyFETXFmv9GUteVLMGVDmxuNJNKoo7Mjcxk5JjO/7TvN0qRkftpxgrTMXAJ9POjfNpxB/4ika4t6eFajBOIoVpOJ1joDeAVAKXUXxr4l0mciRE1wdJPRqb5rsbHkSfvboesEqOe0BgaXkmMyk7j/DEuTklm+4wTnL+UQ4ONB3yvCGRQbQbeWoXh5SAIpqLTRXKnA9VrrzcCnGH0oQojqSmvY9wusm2XsK+ITBD2eMFb3DWjg7OicLtdk5veDqSxJOs7y7Sc4m5GDv7cHfa5owMB2EfSIDsXbQ/ant6a0Zq46gLfl69EYqwYfdHhEQgj7MuXA9m+NJHJqBwREQt+XoONd4F2x1WZrCpNZ8/tBSw1k+wnOXMzGz8ud69s0YFBsBD2jw/DxlARii9KSyWHgfqWUN8ZorquUUj7WCmut19g7OCFEJWRdgD8/N0Zanf8bwlrDTe9B2+Hg4eXs6BzmVFomD3/5J3PuuIr6AcU/skxmzcZDqSzdlsyybSc4fSELX093rmtTn0GxEfSKqS8JpAJKSyYzgQ+AuzCGBb9rpVzesGH57gvhCi6kwB9zjUfmOWjSFQa8Bq36glvNb+efvWIvGw6lMnvFPl68qS0AZrNm85GzLElKZtm2ZE6lZ+Hj6ca1reszsF0kvVuH4edVro1nRRGldcB/rJT6AYjGmLQ4AdhVVYEJIcop9SAkzoE/50NuJrQeBN0ehcadnB1ZlTmVlsmCTUfRGhZu/JteMWEk7DvDsm3JnEjLxMvDjd4xYQyMjeS61vWp4y0JxF5K/U5qrZOBZKXUp8BSrbX0mQjhao7/Cetmw87vQLnDP24zRmaFRTs7sir35i97MJmNuXOZuWbu+3QjXu5uxMeE8a/Y1lzXpgH+kkAcwtaFHsdYO6eUigfu0lrfY7eohKitrO2n3vXhwrPJtYYDK42FFw+uBu9A6PoIXPMgBEY4L/4qlJ1r5q8T6Ww7dp5tx87z55Gz7D6RXqiMp7ti+WM9aB4mKxc7WoVStFKqJcYIrzuBpkAGIMlEiMooaT/17HTY/Bls/dJY56p5b6MGsm4WnEgC/3C4/jmIG2MM9a2hiiaO7cfO89eJdLJNZgACfTzw83LHTYG5yKIeH687lN93IhzH6nIqxQoqFYSxvMpdQGfL4a0YnfRf2roMvVKqPzALo8P+P1rrmUXOPwHcB+RizG25R2t9uNiFCpDlVES1l3rAWIG3pP3U87h7gl8opCdDvVbQbYKxBpaHt/XXVEO2JI52jYJo2zCIdg2DiG0YjLeHouerq8jKNRe7no+HG2ue6l3iyK7arqqWU0Ep5Qb0x0ggNwI+GGt1vQOMBx4rz5BgpZS75bV9gKPABqXU91rrnQWK/QnEaa0zlFIPYszCr70bIojawZb91E05YDbBrf+FmAE1YmRWdq6ZPSfTSTpaeuIY0z0qP3E0rutbbBOpKYu2Ybbyh7FJ60Iju4RjlDYD/nXgDqA+kAkswpgJ/wsQCDxcgft1AvZprQ9Y7vEVMATITyZa65UFyq8HRlXgPkJUL7bspw6QewnaDHJ8PA5Q3sTRrmEQTer6lbn7YN4IrhxTyckkx6RZuPFvJlzXUmonDlRazeRxjPkjy4C7tdZn8k4opSq6HH1DjP1R8hwFriml/L3ADyWdUEqNBcYCNGnSpILhCOECci7VuP3UHZU4SjJ7xV6rtZI8UjtxvNKSyUfALcBA4C9LLeKzqlrsUSk1CogD4ks6r7WeC8wFo8+kKmISwi5yMuHYRji4Fg79Bkc3UJ33U89LHNuOnSfpqGMTR0k2HzlntVaSJ8ek2Xz4rF3uJ0pW2qTF+5VSjwBDMfpMHgAeVErtwWjyqsgH+DGgcYHnjSzHClFKXQ88A8Rrra3toSJE9ZCbBUc3Gonj0FpjL3VTFig3CI+FTvdDyh448KvL76duS+Jo29BxiaMkyx7t4bBrC9uVZzRXBMZQ4NHAFZbD6zGWWVmotc604RoewB7gOowksgG4Q2u9o0CZq4CFQH+t9V5bYpPRXMKl5GbDsU1G4shLHrmZgILwdtCsJ0R1hyZdwDfYeI0to7k8/YxdDato98KCiWPbsfNsO1py4mjXKKjKEoewL3uO5rI5mRQJIA6jtnIbUA84r7UOsfG1A4C3MIYGf6y1fkkp9TywUWv9vVLqF6Adxs6OAEe01oNLu6YkE+FUudlwfLMlefwGR343OspREN7W2LUwqjs07Qq+pfyalDTPBArvp96qj0Pegs2Jo0DykMRR/Tk9mRQIxBMYBIzWWg+1R0AVIclEVClTjrGEyaG1Rr/H379frlE0aGskjqgeRvLwq1u+a6cecPh+6mUljgAfj/yahiSOms1lkomrkGQiHMqUC8lb4OAaS81jPeRcNM7Vv6JAzaMb1Knn3FiLKJo4th87z+7kkhNH24ZBxDaSxFGbVNmkRSFqJVMunNh6ebTVkcTLQ3LD2kD7Oyy1j+5QJ9S5sRZga+IY0y1KEoewO0kmQphNkLz18mirw4nGmlgAoTFGM1OzHtC0O/iHOTdWi/ImjnYNg2haTxKHcBxJJqL2MZvgxLYCySMBsixLy9VrBbG3WJqtuld6b/Sydv2zhSQOUR1IMhE1n9kMJ7dfHm11eB1knjfO1WsJbW++3O8REG7XW5e0619pJHGI6kqSiXA9tu7pYY3ZDKd2GInj4FpL8jhnnKvbHK4YAlE9IaobBEY67G0U3fWv6NpQkjhETSLJRLgWW/b0KDrXwmyGlF2WDnNL8rhkWTojJAra3GipeXSDoEZV9lYKrhll0prp3++gR6swq4mjbWThxNGkrh9ubpI4RPUgQ4OF67B1Fvi4dcaM8vw+j3WQYVmHNLjp5SarqO4Q3Nj6tewsx2TmZFomx89lsis5jReW7CS36E5NXE4csQX25JDEIZxBhgaLmsmWPT1yLsG7nY21rQCCGkOrfpbRVt0gpKlDQtNak3Ypl2PnLnH83CWOn79k+TrTeH7uEifTMovt8leQu5tiULsI3ry1vSQOUeNIMhGuw6Y9PTRoMwx5x6h5hETZ5dY5JjMnzmfmJ4rj5zLzE8exs8a/F7NNhV7j5e5GRLAPkUG+dG0RSsNgHyKDffHzcuefC5Lym7DymMyaH3ec4PTFLNlXQ9Q4kkxE1TPlwvkjcGY/nNl3+ZE3t6MsOheusn3PNK015y/lFKtJ5NcyzmVyMj2Toi2+det4ERnsQ7PQOnRrGUrDYF8ig31pGOJLZLAPoXW8S6xhTFm0DW1lUW3ZV0PUVJJMhGNoDeknCieLM/shdT+kHixcA/EOguDnroUAAA6cSURBVNCW4OZR+hLseYrs6ZGda9QqLieHvGaoy4kjo4RaRaSlJtG9VaiRJCzPI4N9iQzyxdfLvdxvW3b9E7WVJBNROZfOWmoY+4snjrz1qwA8fKBuCwiLgdYDjfkdeQ+/eqAULHnCGLVVSlOXSXnwZ1BfPvlic34TVMqFrGK1ilB/LyKDfWkZ5k/PVmFEBvvk1ywig32pV8fLIf0WsuufqK0kmYiy5VwyRloVTBR5X2ecuVxOuRmjqeq1NDrD67W4nDACG4KbW7FLZ+WaOJGawbFzlzgfOJzr+S+eWE8mWWZ3njreA3NGGpHBPsRHh+U3PeUli4ggH3w8y1+rsAfZ9U/UVjI0WBhMuXDucMk1jLSjhcsGRFiSRIFkUbeF0Rnu4ZVfTGtN6sXsQp3ZRZugUtILb6TZy20L73nNwoNcPLncNGVWHuDuSfrgjwlsd4NM3BPCDmRosKgYrSE9uYQaxn44e7Bwf4VPkLFOVVR3S9JocTlpeBt9Fpk5pvwRUMcOX+L41sOXk8VZo4M7K7fwiCYfTzdL/4QvrWPqW5qdLjdBhQf1xyf99mJ7erhZ9vQIqqJdBoUQ5SM1k5ooI7VwDSM17+sDRfoxfI3lRQrWMCwP7RvCmYycAiOfLndm5z0/fSGr2K3DArwvd2YH+RZrggrx85RahRAuQmomArIvWvox9hfvx7iUermccjcm8tVracwMtySOzKBmHDeFcPx89uVhsvsucXzTJY6f28bxEmoVvp7u+SOg2kQEXh75ZKlZhAf54O3hnL4KIYRzSTJxZaYcOHekSB9GXj/GscJlAyKNRHHFEMx1W5BepynJ7g05ZA7laJrJmF+RYiSMY2cvcebi9kIvVwrqW2oVV0QG0ueKBkQGXR4q2zDYl2CpVQghrJBk4mxmc8n9GKn74eyhIv0YwfD/7Z17jFTVHcc/X3bZB499IOwuuwiCShGMWp/VoNE2BrUxaIvRWoNFW5M2NphWo63xEZs2PtKmDx8tpa3aGo3xFdtqrI+qrVoVLQgoyMIC8pLiorDCsuzur3+cM+zdcdHdmdmZWfP7JDdz5txzz/3umTP7m3N+5/7O2EPpnDiTj0ZOZOvwCaynkXc761i3U2EKasVuNn3UTkdnN7AlHjCirGSfYZjRWEVjderhu2Ao6qsqKCv95Gorx3Gc/uDGJF/sav3k6CL1EF8isKGVVtJZM5m20VP5YOxX2DCskeauBpbvGceqtjI2bW6ntbkjUfFepI3Uj66gsaaCw5uqmTWj4RNTUNWVPqpwHGfwcGOSSzo+7jEQn/Bj9DxXYCph96gD2V4xkS1jZ9NijazYW8fiXQewdMcI9rT1rnZkWQlNtaU01pRzxISa6MzucXA3VFcwvMRHFY7jFA43JgOlay9sX9e3H2Pnpl5F28rr2Vo2gfWlM1lVWcfS9nEsbR/HezaOzt2h6YcJ6quib2JiJcfEkUTyae2qilIfVTiOU9R8PozJpsV0/6wpPIvQ3934Po3u7mAY0vwY3dua0YfrkPU8TNdWUs2mYeNZ3T2Vd7pm0tzVQIuNZ63Vs7u9glHlpcEw1AWDMaeXoaigvspHFY7jDH0+H8+ZNJbYostGwbDhUDK879340jHr5cfo3tZMx9ZVdG9bRfmOtZR0te8r2k45a2nYZyhausPreo2nomrcvuWyPSufKhKjiuGD/Nc7juNkxpB+zkTSGcCvgBJgoZndnHa+HLgXOAb4ADjfzNb2q/LuuNXrg3Phuy+FEcqeNmhdTfuWd2nbvJLOraso3b6akW3rqOzase/SLitho9XRYg202JdpsfFsGd5ER9VkymubGF87gsaaSibWVHJiNBR1o8sp9VGF4zhOfo2JpBLgDuB0YAPwuqTHzeztRLFLge1mdoikC4BbgPMHcp/uvbvZeefpdJtR2xUCEVbEY6MdQHN3A2s5ntbyA/l41EF0jTmYirGTaRgzmqaaSk6pqeT8mgofVTiO4/STfI9MjgeazWwNgKQHgNlA0pjMBm6M6YeA2yXJBjAfNwxjVGcrz5Wdxs6qSXRUT2HYuEMY2TCVhrG1HFxTyYmjKyjxrVMdx3FyQr6NSRPwXuL9BuCE/ZUxs05JHwEHANuShSRdBlwGcEAlHLsgbT0t8MbmR9/ImfKBM5Y0zUVCMepyTf3DNfWfYtRVjJq+kKuKhuxqLjNbACwAkLRo266unDiRcoWkRblybOWSYtTlmvqHa+o/xairWDXlqq58e483Agcm3k+IeX2WkVQKVBMc8Y7jOE6Rkm9j8jpwqKTJksqAC4DH08o8Dlwc03OA5wbiL3Ecx3HyT16nuaIP5HLgKcLS4D+a2XJJNwGLzOxx4A/AnyU1A60Eg/NZLBg00ZlTjJqgOHW5pv7hmvpPMer6XGv6XDy06DiO4xQWf+LOcRzHyRo3Jo7jOE7WFKUxkXSgpH9KelvScknzY/4YSU9LWhVfa2P+NEmvSNoj6cq0us6QtFJSs6Rr8qzrm5LekrRU0suSjsy1rgw0zY6aFktaJGlmoq6LY/lVki7e3z1zrSlx3XGSOiXNKbQmSadK+ii202JJ1yfqKshnl9C1OJZ/IdeaMtEl6apEOy2T1CVpTC51ZaCpWtJfJS2J5ecl6ipUn6qV9Gj8/r0m6fBEXYPdTufF992Sjk275kfxvislzcpYk5kV3QGMB46O6dHAu8B04Fbgmph/DXBLTNcBxwE/Ba5M1FMCrAamAGXAEmB6HnWdBNTG9JnAq7nWlYGmUfT4yo4AVsT0GGBNfK2N6dp8aEq0yXPAE8CcQmsCTgX+1kc9hfzsagjRIiam+n0x9PO0a88mrMAsdFv9OJEeR1jMU1bgPnUbcENMTwOezWM7HUZ4QPF54NhE+enxfuXA5KijJBNNRTkyMbPNZvZmTO8E3iE8GT8buCcWuwc4J5bZamavA3vTqtoXvsXMOoBU+JZ86XrZzFK7Yv2H8FxNTnVloKnNYi8CRgKp9CzgaTNrjZqfBs7Ih6bI94GHga2JvEJr6ouCfXbAhcAjZrY+XpNqq4L28zS+Adyfa10ZaDJgtCQRfkC1Ap0Utk9NJ/xgwsxWAAdJqicP7WRm75jZyj4umQ08YGZ7zKwFaI56BqypKI1JEkkHAV8EXgXqzWxzPLUFqP+My/sK39JUIF2XAk8Opq7+apJ0rqQVwN+BSwqtSVITcC5wV9rlBW0n4MQ4TfKkpBlFoGkqUCvpeUlvSJo7mJoGoCtVdgThH/PDg6mrn5puJ/wa3wQsBeabWXeBNS0BvhbLHw9MIvzAzIem/bG/ew9YU1EbE0mjCB3zCjPbkTwXf10XZF3zQHVJOo1gTK4uBk1m9qiZTSP8YvpJEWj6JXB1/LIPKgPQ9CYwycyOBH4DPFYEmkoJWzN8lfAL+zpJU4tAV4qzgZfMrLUINM0CFgONwFGEgLFVBdZ0M1AjaTFhJP5foItB4NM0DRZFa0wkDSc0xn1m9kjMfl/S+Hh+PL2nRPqiP+FbBlWXpCOAhcBsM0uFhcmprkzbysxeBKZIGltgTccCD0haS4h6cKekcwqpycx2mFlbTD8BDC+CdtoAPGVmH5vZNuBF4Mhca8pAV4oL6JniIte6BqhpHmFK0MysGWgh+CkK3afmmdlRwFyCL2dNnjTtj/3de+CaLEMn3WAegAgbZP0yLf82eju2bk07fyO9HfClhA9rMj1OpBn50gVMJMxBnpRWPme6MtB0CD0O+KNjBxHBIdlCcErWxvSYfH5+Mf9uejvgC6IJaEi00/HA+lhHIT+7w4Bno4YRwDLg8EL38/i+muCXGFkk/fwu4MaYro/9fGyB+1QNUBbT3wHuzVc7Jc4/T28H/Ax6O+DXEJzvA9aUUWcb7AOYSRgavkUYqi4GziKEon8WWAU8k+oEhC/+BmAH8GFMV8VzZxFWNKwGrs2zroXA9kTZRYm6cqIrA01XA8tjuVeAmYm6LiEYv2ZgXr40pV17N9GYFFITcHlspyWExRMnJeoqyGcXr7mKsKJrGWEKI6eastD1LYIjN72uQvXzRuAfBH/JMuCiIuhTJ8a2WAk8QmIVWR7a6VzC/8U9wPuEEW7qmmvjfVcCZ2aqycOpOI7jOFlTtD4Tx3EcZ+jgxsRxHMfJGjcmjuM4Tta4MXEcx3Gyxo2J4ziOkzVuTBynDyQ9JKk1xk5KP3dqjL46vxDaHKcY8aXBjtMH0Yi8TYiAe14iv5Kwhn8rcLLlIQSM4wwFfGTiOH1gZu8D84E5MbRLihsJoSUuyZchkVQiqTQf93KcTHFj4jj7wcz+QoisfEfcbOlo4IeEMB29wnlLOj9ueLRL0nZJD8RoyMkycyW9IOl/knbGyL8XppWpkGSSrpd0naR1QAdw6OD+tY6THT7N5TifQjQIy4FHCdFnO4EvmVlXoswVwC+A3xOiC9cANxHCWhxlZrtiuRsI02PN8dLTCOFtLjWzu2OZCmA3IXT6SkLE4nbgNesJFOo4RYcbE8f5DCR9m2Ao9gLHmNnSxLkaQhDBe8zse4n8qQSfy+Vm9ts+6hxGmBn4EzDVzE6I+Sljsi7mdwzaH+Y4OcSnuRznMzCzhcBm4LGkIYmcTIjge5+k0tRBiLi6BjglVVDSNEkPStpEGOHsBS4ibKeazhNuSJyhhDv1HKd/dMQjnbr4+u/9XNcC+0YwzxDCtF8V8zuAKwh7uKSzuY88xyla3Jg4Tnak/BgXEkKOp5Pa5e5kwran55jZotTJuJFRX/j8szOkcGPiONnxIsHHMcXM7v+UciPi695UhqQ6wp4RjjPkcWPiOFlgZq2SrgF+LqkReArYSRiFnAY8aWYPAf8CPgZ+J+kmoAq4nrBR0YSCiHecHOLGxHGyxMx+HZ8H+QFhb+8SwgqvFwg7/WFmmyR9HbiVsD/3BsJy4kkEv4njDGl8abDjOI6TNb402HEcx8kaNyaO4zhO1rgxcRzHcbLGjYnjOI6TNW5MHMdxnKxxY+I4juNkjRsTx3EcJ2vcmDiO4zhZ83/2hXAPoQNWqgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(6,4))\n",
    "\n",
    "plt.plot([2010] + medium_years, [0] + sum_medium_results, label=\"Medium SLR scenario\", marker=\"^\", markersize=10)\n",
    "plt.plot([2010] + high_years, [0] + sum_high_results, label=\"High SLR scenario\", marker=\"o\", markersize=10)\n",
    "\n",
    "med_line = matplotlib.lines.Line2D([], [], color='#1f77b4', marker='^', linestyle='-', markersize=8, label='Medium SLR scenario')\n",
    "high_line = matplotlib.lines.Line2D([], [], color='#ff7f0e', marker='o', linestyle='-', markersize=8, label='High SLR scenario')\n",
    "\n",
    "plt.legend(loc=\"best\", fontsize=16, handles=[high_line, med_line])\n",
    "\n",
    "plt.xlim([2010,2100])\n",
    "plt.ylim([0, 1.3e7])\n",
    "\n",
    "plt.xlabel(\"Year\", fontsize=16)\n",
    "plt.ylabel(\"Affected Population\", fontsize=16)\n",
    "\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save results per county"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = open(\"data/processed/affected_population_medium.csv\",\"w\")\n",
    "f.write(\"County FIPS,\")\n",
    "for i,year in enumerate(medium_years):\n",
    "    f.write(\"Total Population %d,Affected Population %d,\" % (year, year))\n",
    "f.write(\"\\n\")\n",
    "for county_idx in range(len(migration_counties)):\n",
    "    f.write(\"%s,\" % (migration_counties[county_idx]))\n",
    "    for i,year in enumerate(medium_years):\n",
    "        f.write(\"%f,\" % (medium_results_county[i][0][county_idx]))\n",
    "        f.write(\"%f,\" % (medium_results_county[i][1][county_idx]))\n",
    "    f.write(\"\\n\")\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = open(\"data/processed/affected_population_high.csv\",\"w\")\n",
    "f.write(\"County FIPS,\")\n",
    "for i,year in enumerate(high_years):\n",
    "    f.write(\"Total Population %d,Affected Population %d,\" % (year, year))\n",
    "f.write(\"\\n\")\n",
    "for county_idx in range(len(migration_counties)):\n",
    "    f.write(\"%s,\" % (migration_counties[county_idx]))\n",
    "    for i,year in enumerate(high_years):\n",
    "        f.write(\"%f,\" % (high_results_county[i][0][county_idx]))\n",
    "        f.write(\"%f,\" % (high_results_county[i][1][county_idx]))\n",
    "    f.write(\"\\n\")\n",
    "f.close()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}