Compare the area burnt by Australia's fires to the area of countries and plot out a quick comparison

Area of Australian fires vs countries.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Area of Australia's 2019/2020 bushfires vs a country"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'retina'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import os.path\n",
    "\n",
    "from descartes.patch import PolygonPatch\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import optimize\n",
    "from shapely import geometry, ops"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data prep\n",
    "\n",
    "First grab the data. [Natural Earth](http://naturalearthdata.com) low resolution is perfect for boundaries to startm and the [CIA world factbook](https://www.cia.gov/library/publications/resources/the-world-factbook/) is a handy source for a list of country areas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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>featurecla</th>\n",
       "      <th>scalerank</th>\n",
       "      <th>LABELRANK</th>\n",
       "      <th>SOVEREIGNT</th>\n",
       "      <th>SOV_A3</th>\n",
       "      <th>ADM0_DIF</th>\n",
       "      <th>LEVEL</th>\n",
       "      <th>TYPE</th>\n",
       "      <th>ADM0_A3</th>\n",
       "      <th>GEOU_DIF</th>\n",
       "      <th>...</th>\n",
       "      <th>NAME_KO</th>\n",
       "      <th>NAME_NL</th>\n",
       "      <th>NAME_PL</th>\n",
       "      <th>NAME_PT</th>\n",
       "      <th>NAME_RU</th>\n",
       "      <th>NAME_SV</th>\n",
       "      <th>NAME_TR</th>\n",
       "      <th>NAME_VI</th>\n",
       "      <th>NAME_ZH</th>\n",
       "      <th>geometry</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ADMIN</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>New Zealand</th>\n",
       "      <td>Admin-0 country</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>New Zealand</td>\n",
       "      <td>NZ1</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>Country</td>\n",
       "      <td>NZL</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>뉴질랜드</td>\n",
       "      <td>Nieuw-Zeeland</td>\n",
       "      <td>Nowa Zelandia</td>\n",
       "      <td>Nova Zelândia</td>\n",
       "      <td>Новая Зеландия</td>\n",
       "      <td>Nya Zeeland</td>\n",
       "      <td>Yeni Zelanda</td>\n",
       "      <td>New Zealand</td>\n",
       "      <td>新西兰</td>\n",
       "      <td>MULTIPOLYGON (((173.11533 -41.27930, 173.23086...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Nauru</th>\n",
       "      <td>Admin-0 country</td>\n",
       "      <td>3</td>\n",
       "      <td>6</td>\n",
       "      <td>Nauru</td>\n",
       "      <td>NRU</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>Sovereign country</td>\n",
       "      <td>NRU</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>나우루</td>\n",
       "      <td>Nauru</td>\n",
       "      <td>Nauru</td>\n",
       "      <td>Nauru</td>\n",
       "      <td>Науру</td>\n",
       "      <td>Nauru</td>\n",
       "      <td>Nauru</td>\n",
       "      <td>Nauru</td>\n",
       "      <td>諾魯</td>\n",
       "      <td>POLYGON ((166.95840 -0.51660, 166.93896 -0.550...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>2 rows × 94 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                  featurecla  scalerank  LABELRANK   SOVEREIGNT SOV_A3  \\\n",
       "ADMIN                                                                    \n",
       "New Zealand  Admin-0 country          1          2  New Zealand    NZ1   \n",
       "Nauru        Admin-0 country          3          6        Nauru    NRU   \n",
       "\n",
       "             ADM0_DIF  LEVEL               TYPE ADM0_A3  GEOU_DIF  ...  \\\n",
       "ADMIN                                                              ...   \n",
       "New Zealand         1      2            Country     NZL         0  ...   \n",
       "Nauru               0      2  Sovereign country     NRU         0  ...   \n",
       "\n",
       "            NAME_KO        NAME_NL        NAME_PL        NAME_PT  \\\n",
       "ADMIN                                                              \n",
       "New Zealand    뉴질랜드  Nieuw-Zeeland  Nowa Zelandia  Nova Zelândia   \n",
       "Nauru           나우루          Nauru          Nauru          Nauru   \n",
       "\n",
       "                    NAME_RU      NAME_SV       NAME_TR      NAME_VI NAME_ZH  \\\n",
       "ADMIN                                                                         \n",
       "New Zealand  Новая Зеландия  Nya Zeeland  Yeni Zelanda  New Zealand     新西兰   \n",
       "Nauru                 Науру        Nauru         Nauru        Nauru      諾魯   \n",
       "\n",
       "                                                      geometry  \n",
       "ADMIN                                                           \n",
       "New Zealand  MULTIPOLYGON (((173.11533 -41.27930, 173.23086...  \n",
       "Nauru        POLYGON ((166.95840 -0.51660, 166.93896 -0.550...  \n",
       "\n",
       "[2 rows x 94 columns]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "world = gpd.read_file(\n",
    "    os.path.expanduser('~/Developer/spatial_data/naturalearth/natural_earth_vector.gpkg'),\n",
    "    layer='ne_50m_admin_0_countries'\n",
    ")\n",
    "world.set_index('ADMIN', inplace=True)\n",
    "world.sample(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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>area</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>name</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Tonga</th>\n",
       "      <td>747.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Bahamas, The</th>\n",
       "      <td>13880.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                 area\n",
       "name                 \n",
       "Tonga           747.0\n",
       "Bahamas, The  13880.0"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "country_areas = pd.read_fwf(\n",
    "    'https://www.cia.gov/library/publications/resources/the-world-factbook/fields/rawdata_279.txt',\n",
    "    widths=[8, 57, 21],\n",
    "    header=None,\n",
    "    skiprows=2,\n",
    "    names=['id', 'name', 'area'],\n",
    "    thousands=',',\n",
    "    index_col='id'\n",
    ")\n",
    "country_areas.set_index('name', inplace=True)\n",
    "country_areas.sample(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note** As this is a thought exercise, not going to worry about missing matches - it's an all right starting point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "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>featurecla</th>\n",
       "      <th>scalerank</th>\n",
       "      <th>LABELRANK</th>\n",
       "      <th>SOVEREIGNT</th>\n",
       "      <th>SOV_A3</th>\n",
       "      <th>ADM0_DIF</th>\n",
       "      <th>LEVEL</th>\n",
       "      <th>TYPE</th>\n",
       "      <th>ADM0_A3</th>\n",
       "      <th>GEOU_DIF</th>\n",
       "      <th>...</th>\n",
       "      <th>NAME_NL</th>\n",
       "      <th>NAME_PL</th>\n",
       "      <th>NAME_PT</th>\n",
       "      <th>NAME_RU</th>\n",
       "      <th>NAME_SV</th>\n",
       "      <th>NAME_TR</th>\n",
       "      <th>NAME_VI</th>\n",
       "      <th>NAME_ZH</th>\n",
       "      <th>geometry</th>\n",
       "      <th>area</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Cook Islands</th>\n",
       "      <td>Admin-0 country</td>\n",
       "      <td>3</td>\n",
       "      <td>4</td>\n",
       "      <td>New Zealand</td>\n",
       "      <td>NZ1</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>Dependency</td>\n",
       "      <td>COK</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>Cookeilanden</td>\n",
       "      <td>Wyspy Cooka</td>\n",
       "      <td>Ilhas Cook</td>\n",
       "      <td>Острова Кука</td>\n",
       "      <td>Cooköarna</td>\n",
       "      <td>Cook Adaları</td>\n",
       "      <td>Quần đảo Cook</td>\n",
       "      <td>库克群岛</td>\n",
       "      <td>POLYGON ((-159.74053 -21.24922, -159.77256 -21...</td>\n",
       "      <td>236.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>El Salvador</th>\n",
       "      <td>Admin-0 country</td>\n",
       "      <td>1</td>\n",
       "      <td>6</td>\n",
       "      <td>El Salvador</td>\n",
       "      <td>SLV</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>Sovereign country</td>\n",
       "      <td>SLV</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>El Salvador</td>\n",
       "      <td>Salwador</td>\n",
       "      <td>El Salvador</td>\n",
       "      <td>Сальвадор</td>\n",
       "      <td>El Salvador</td>\n",
       "      <td>El Salvador</td>\n",
       "      <td>El Salvador</td>\n",
       "      <td>萨尔瓦多</td>\n",
       "      <td>POLYGON ((-89.36260 14.41602, -89.33726 14.411...</td>\n",
       "      <td>21041.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>2 rows × 95 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                   featurecla  scalerank  LABELRANK   SOVEREIGNT SOV_A3  \\\n",
       "Cook Islands  Admin-0 country          3          4  New Zealand    NZ1   \n",
       "El Salvador   Admin-0 country          1          6  El Salvador    SLV   \n",
       "\n",
       "              ADM0_DIF  LEVEL               TYPE ADM0_A3  GEOU_DIF  ...  \\\n",
       "Cook Islands         1      2         Dependency     COK         0  ...   \n",
       "El Salvador          0      2  Sovereign country     SLV         0  ...   \n",
       "\n",
       "                   NAME_NL      NAME_PL      NAME_PT       NAME_RU  \\\n",
       "Cook Islands  Cookeilanden  Wyspy Cooka   Ilhas Cook  Острова Кука   \n",
       "El Salvador    El Salvador     Salwador  El Salvador     Сальвадор   \n",
       "\n",
       "                  NAME_SV       NAME_TR        NAME_VI NAME_ZH  \\\n",
       "Cook Islands    Cooköarna  Cook Adaları  Quần đảo Cook    库克群岛   \n",
       "El Salvador   El Salvador   El Salvador    El Salvador    萨尔瓦多   \n",
       "\n",
       "                                                       geometry     area  \n",
       "Cook Islands  POLYGON ((-159.74053 -21.24922, -159.77256 -21...    236.0  \n",
       "El Salvador   POLYGON ((-89.36260 14.41602, -89.33726 14.411...  21041.0  \n",
       "\n",
       "[2 rows x 95 columns]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "world = world.join(country_areas, how='inner')\n",
    "world.sample(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Performing the calculations and producing the partial polygon\n",
    "\n",
    "Ideally we'd like to know  what percentage of the area of the countries have been burned by Australia's fires. If we can it'd be nice to get a polygons that represent this number visually. We can get the total area burned from the [SBS report](https://www.sbs.com.au/news/the-numbers-behind-australia-s-catastrophic-bushfire-season):\n",
    "\n",
    "- NSW/ACT: 3.6 million ha\n",
    "- Western Australia: 1.2 million ha\n",
    "- Victoria: 780,000ha\n",
    "- South Australia: 60,000ha\n",
    "- Queensland: 250,000ha\n",
    "- Tasmania: 8000ha\n",
    "\n",
    "For a total of 5,898,000ha burned (so far), or 58,980km<sup>2</sup> burned. To put that in proportion, in this dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "burned_area = 58_980  # km**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "55% of countries have greater area than the burned area in these fires\n"
     ]
    }
   ],
   "source": [
    "print(\n",
    "    '{:.0f}% of countries have greater area than the burned area in these fires'.format(\n",
    "        (world['area'] > burned_area).mean() * 100\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So let's pick a country. For whatever reason, Belgium is the big comparison of the day:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"100.0\" height=\"100.0\" viewBox=\"2.3713203124999995 49.35730664062499 4.146714843750011 2.287388671875007\" preserveAspectRatio=\"xMinYMin meet\"><g transform=\"matrix(1,0,0,-1,0,101.00200195312499)\"><path fill-rule=\"evenodd\" fill=\"#66cc99\" stroke=\"#555555\" stroke-width=\"0.08293429687500023\" opacity=\"0.6\" d=\"M 4.226171875000006,51.386474609375 L 4.304492187500017,51.3615234375 L 4.373730468750011,51.356005859374996 L 4.404003906250011,51.367089843749994 L 4.384765625,51.427587890625 L 4.44091796875,51.459814453125 L 4.50341796875,51.47470703125 L 4.531640625000023,51.448583984375 L 4.588769531250023,51.421923828124996 L 4.633984375000011,51.421728515625 L 4.755664062500017,51.491113281249994 L 4.7841796875,51.477392578125 L 4.810546875,51.452734375 L 4.816015625000006,51.4328125 L 4.820703125000023,51.412060546875 L 4.848046875000023,51.403271484375 L 4.943945312500006,51.407763671874996 L 4.992578125000023,51.445361328124996 L 5.030957031250011,51.469091796875 L 5.059472656250023,51.453125 L 5.073437500000011,51.4068359375 L 5.099902343750017,51.346484374999996 L 5.214160156250017,51.278955078124994 L 5.310839843750017,51.259716796875 L 5.429785156250006,51.272998046874996 L 5.476855468750017,51.285058593749994 L 5.5087890625,51.275 L 5.540429687500023,51.239306640624996 L 5.608789062500023,51.1984375 L 5.752343750000023,51.169482421874996 L 5.796484375000006,51.153076171875 L 5.8271484375,51.125634765624994 L 5.818261718750023,51.08642578125 L 5.749804687500017,50.98876953125 L 5.740820312500006,50.959912109375 L 5.75,50.950244140624996 L 5.736621093750017,50.93212890625 L 5.647558593750006,50.866650390625 L 5.639453125000017,50.843603515625 L 5.669140625000011,50.805957031249996 L 5.693554687500011,50.774755859375 L 5.693652343750017,50.774658203125 L 5.7470703125,50.759570312499996 L 5.79736328125,50.754541015624994 L 5.892461122495295,50.752556857983514 L 5.993945312500017,50.750439453125 L 6.005957031250006,50.732226562499996 L 6.119433593750017,50.679248046874996 L 6.154492187500011,50.637255859374996 L 6.235937500000006,50.5966796875 L 6.16845703125,50.545361328125 L 6.1787109375,50.522509765624996 L 6.203027343750023,50.49912109375 L 6.294921875,50.485498046874994 L 6.340917968750006,50.4517578125 L 6.343652343750023,50.400244140625 L 6.364453125000011,50.316162109375 L 6.175097656250017,50.232666015625 L 6.121289062500011,50.13935546875 L 6.116503906250017,50.120996093749994 L 6.110058593750011,50.123779296875 L 6.089062500000011,50.154589843749996 L 6.054785156250006,50.154296875 L 5.976269531250011,50.1671875 L 5.866894531250011,50.082812499999996 L 5.8173828125,50.0126953125 L 5.7880859375,49.961230468749996 L 5.744042968750023,49.91962890625 L 5.735253906250023,49.875634765624994 L 5.740820312500006,49.857177734375 L 5.725781250000011,49.833349609375 L 5.725000000000023,49.80830078125 L 5.787988281250023,49.75888671875 L 5.8037109375,49.732177734375 L 5.88037109375,49.644775390625 L 5.856542968750006,49.612841796874996 L 5.837597656250011,49.5783203125 L 5.8154296875,49.55380859375 L 5.789746093750011,49.53828125 L 5.71044921875,49.539208984374994 L 5.610058593750011,49.528222656249994 L 5.542382812500023,49.511035156249996 L 5.50732421875,49.510888671874994 L 5.434667968750006,49.554492187499996 L 5.353515625,49.61982421875 L 5.301953125000011,49.6509765625 L 5.27880859375,49.6779296875 L 5.215039062500011,49.689257812499996 L 5.124121093750006,49.721484374999996 L 5.06103515625,49.75654296875 L 5.006933593750006,49.778369140624996 L 4.930566406250023,49.7892578125 L 4.867578125000023,49.788134765624996 L 4.84912109375,49.847119140625 L 4.841503906250011,49.914501953125 L 4.7900390625,49.9595703125 L 4.860546875000011,50.135888671874994 L 4.818652343750017,50.153173828125 L 4.772851562500023,50.139062499999994 L 4.706640625000006,50.097070312499994 L 4.675097656250017,50.046875 L 4.656152343750023,50.00244140625 L 4.545019531250006,49.960253906249996 L 4.368750000000006,49.944970703124994 L 4.176074218750017,49.960253906249996 L 4.149316406250023,49.971582031249994 L 4.137011718750017,49.98447265625 L 4.136816406250006,50.0 L 4.150292968750023,50.023876953125 L 4.183886718750017,50.05283203125 L 4.192187500000017,50.094140624999994 L 4.15771484375,50.1298828125 L 4.13525390625,50.143798828125 L 4.144140625000006,50.17841796875 L 4.169628906250011,50.22177734375 L 4.174609375000017,50.246484374999994 L 4.044140625000011,50.321337890624996 L 3.94970703125,50.3359375 L 3.8581054687500114,50.33857421875 L 3.78857421875,50.346972656249996 L 3.748046875,50.343505859375 L 3.7188476562500057,50.321679687499994 L 3.6893554687500227,50.306054687499994 L 3.667285156250017,50.3248046875 L 3.626757812500017,50.457324218749996 L 3.5954101562500114,50.477343749999996 L 3.4769531250000227,50.499462890625 L 3.316210937500017,50.507373046874996 L 3.2733398437500227,50.531542968749996 L 3.249804687500017,50.591162109375 L 3.2349609375000057,50.662939453125 L 3.1820312500000227,50.731689453125 L 3.1548828125000057,50.748925781249994 L 3.1068359375000227,50.779443359374994 L 3.0228515625000227,50.766894531249996 L 2.9219726562500057,50.72705078125 L 2.8624023437500057,50.716015625 L 2.8397460937500227,50.711767578125 L 2.7593750000000057,50.750634765624994 L 2.6691406250000114,50.811425781249994 L 2.5967773437500057,50.875927734375 L 2.5792968750000114,50.911767578124994 L 2.6014648437500227,50.9552734375 L 2.5748046875000057,50.988574218749996 L 2.5360351562500227,51.04951171875 L 2.52490234375,51.097119140625 L 2.9601562500000114,51.265429687499996 L 3.2251953125000057,51.351611328124996 L 3.35009765625,51.377685546875 L 3.3800781250000114,51.29111328125 L 3.40283203125,51.263623046875 L 3.4325195312500227,51.245751953124994 L 3.471972656250017,51.242236328124996 L 3.51708984375,51.263623046875 L 3.5802734375000114,51.286181640624996 L 3.6818359375000114,51.275683593749996 L 3.755664062500017,51.254833984375 L 3.7819335937500114,51.233203124999996 L 3.8307617187500114,51.21259765625 L 3.9020507812500114,51.207666015624994 L 4.0400390625,51.2470703125 L 4.172558593750011,51.307080078125 L 4.21142578125,51.34873046875 L 4.226171875000006,51.386474609375 z\" /></g></svg>"
      ],
      "text/plain": [
       "<shapely.geometry.polygon.Polygon at 0x1024826490>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "world.loc['Belgium', 'geometry']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A nice simple shape. So how many Belgiums?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.931996855345912"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "burned_area / world.loc['Belgium', 'area']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we want to be able to plot Belgium 1 time, and partially plot the remainder. How can we get a partial fill? Well, we have the proportion, and we can use some simple optimisation to work out the filled fraction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def clip_at_horizontal_line(geom, y_offset):\n",
    "    \"\"\"\n",
    "    Returns the geometry clipped to the horizontal line defined by the y offset.\n",
    "    \n",
    "    First, calculate the bounding box, add the y_offset to the minimum y value to define the new rectangle, then\n",
    "    intersect that rectangle with the original piece of data. Returns the clipped geometry.\n",
    "    \"\"\"\n",
    "    minx, miny, maxx, maxy = geom.bounds\n",
    "    \n",
    "    clipper = geometry.Polygon([\n",
    "        [minx, miny],\n",
    "        [minx, miny + y_offset],\n",
    "        [maxx, miny + y_offset],\n",
    "        [maxx, miny]\n",
    "    ])\n",
    "    \n",
    "    return clipper.intersection(geom)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def minimise_area_proportional_difference(y_offset, proportion, geom):\n",
    "    \"\"\"\n",
    "    Calculate the absolute value of the proportional difference in areas of the clipped geom and the original.\n",
    "    \"\"\"\n",
    "    clipped_geom = clip_at_horizontal_line(geom, y_offset)\n",
    "    \n",
    "    found_proportion = clipped_geom.area / geom.area\n",
    "    \n",
    "    return abs(found_proportion - proportion)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9319968553459119"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "target_proportion = burned_area / world.loc['Belgium', 'area'] % 1\n",
    "target_proportion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "_, miny, _, maxy = world.loc['Belgium', 'geometry'].bounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "solution = optimize.minimize_scalar(\n",
    "    minimise_area_proportional_difference,\n",
    "    method='Bounded',\n",
    "    bounds=(0, maxy - miny),\n",
    "    args=(target_proportion, world.loc['Belgium', 'geometry'])\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.7179691221651343"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solution.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "clipped_belgium = clip_at_horizontal_line(world.loc['Belgium', 'geometry'], solution.x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting out some countries\n",
    "\n",
    "Let's start with Belgium:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def set_limits(ax, x0, y0, xN, yN):\n",
    "    ax.set_xlim(x0, xN)\n",
    "    ax.set_ylim(y0, yN)\n",
    "    ax.set_aspect(\"equal\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABHgAAAEiCAYAAACC32V1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdeZycVZX/8c/t7qSzr52kk04IhIRAQsKOICCKCuKCIjgzijuos7nMzDPjNlqW66gPzm8cHZ1xV3RmFEUhIKMCIntCIBtLAoGEJJ10Qvalk/Ryf3+cqu7qWrqru6ue2r7v16te3XWfWk4gnb517rnnOu89IiIiIiIiIiJSuepKHYCIiIiIiIiIiAyPEjwiIiIiIiIiIhVOCR4RERERERERkQqnBI+IiIiIiIiISIVTgkdEREREREREpMIpwSMiIiIiIiIiUuGU4BERERERERERqXBK8IiIiIiIiIiIVDgleEREREREREREKpwSPCIiIiIiIiIiFU4JHhERERERERGRCqcEj4iIiIiIiIhIhVOCR0RERERERESkwjWUOgARqUyhcw44Ebgk5eaAywPvN0fw/guAELhqgId6IAi8/1qxY0oXOjcSeBXwu8D7zqjfX0RERKRY4vH4q4B/BU5Pu7QB+Bfgplgs1hF5YCI1zHnvSx2DiFSA0Lk6YDF9EzotWR66FUvyPFWkOCYBnwI+CIwYxFM/D3w6iOgfvdC5S4BvA4uAJ4EAuDOq9xcREZHqFI/H64DT6J2PvQT4WCwWuzmi9z8FW2R7A8CYMWNYsmQJx48fZ/369Rw5ciT50N8DV8discNRxJUW4whgRiwW2xr1e4uUkhI8ItKv0LnpwGeBPwMm5/m03cBrA++XFzCOBuB9wOeAqdkeswjLOD0K7M3+Mt8APhx4312ouNKFzjUBXwHek+Xy74F/CLxfW6z3FxERkeqSSFacTW9C52JgStrDuoG/jMVi3yliHJPpXWRraGho4NRTT2X+/PnU19dbEN3dbNmyhbVr13Ls2DGAh4DXxWKxHFOzosR5MbbItgD4d+DzsVhsX1TvL1JKSvCISFaJ7UUfwn6RTxjCSxwG3hR4/4cCxHI58DWsgijDLOAVwMzE/W7gaSybkqUu+CfAewu9ZSpR4fRu4KtkTrpSdQPfAz6gah4RERHJJR6Pz8AWtq4DxqReGzVqFE1NTTQ1NdHe3s769euTlz4Wi8W+XOA4GoD3Ywt+UwFOPPFEFi9ezKhRo7I+5+DBg9x///3Jap41wBWxWGxHIePKEudU4MvA9WmXdgNx4NvaMibVTgkeEekj0VvnDcCNwPz+HtuAJVVagIPAE9kfdhvwL4H3Dw4hloWJOF6X7fp44FJgIdb8J10r8CvgaOalNwTeLxtsPP3EORv4KfCyQTzt3MD7lYWKQURERKpDPB5vpHeRbTzAuHHjehI6U6dOZezYsdiUzWzcuJFVq1Yl794IfDIWix0rQCyXY312FgE0NTVxxhlnMGnSpAGfe+TIEe677z4OHToE8CywJBaLZZmWDTtGB7wT2zbW5Jxj4cKFNDc3s27dOl588cXkQx8HzonFYvoALFVLCR4R6RE6dzr2S/xV2a43Ysmc2YnbDKA+cc0D9wOP5H75+7BVlTsGqlwJnZsCfBr4G7I0g2/ANpufy8BNeHYBN2PlRCn+PfD+QwM8NS+hc68C/htoSr82NRHnWmBL5lPjgfefKUQMIiIiUvkSiYqrsATNyQDNzc0sWbKECRMGLqbesmULK1asIDHN2o7N6f4zFosdGEIsp2IJk9cBjB07liVLljBr1qw+iaWBHDt2jLvvvjtZyfOqWCx212BjGSDOacBNwOUA06ZN48wzz+z57+W9Z/v27Tz66KN0dHQAnBmLxVYXMgaRcqJTtEQk2Tfms8AHgLr0643AhcBZ9CZ00jlsU/go4N7sD0nuG18bOvdl4H/Tt0mFzo1IxBAnxzanZJfncf3/kXpMw7JVv+k7fEWeT88psSXr41jpdJ+ZTgP23+tc7L/XacByLMOV4irgM8ONQ0RERCpfPB5fgiVkXgkwfvx4li5dSnNzc96vMWfOHEaOHMnatWvZv3//TKwn4Cfj8fi3gH/LZ4tUPB7vs8iWrc/OYDQ2NjJ79mw2bNgAloQpWIInHo9fAPwCmD1y5EiWLl3KCSec0CcB5Zxj1qxZtLS0sGnTJrAqdSV4pGqpgkekxoXOnQQ8BmTU2jpgKXARaRu/B7AReBBo6/9he7GdXT7lNg7LyWRowfrs5D/N6XUM+CbW/CbFnMD7IZ2sEDo3Gevlk7F1bC7wajL/Yx4GvpX5UkOOQURERCpfIqHyeRKLbCNGjGDRokXMmzePurqMNbe8eO9pa2tj/fr1qduTjgE/BMJYLPZsljgyFtlOOukkFi1alLPPTr7a2tq4//77AR6PxWJnD+vF6Kl0+lus0mnE1KlTeclLXsLo0aNzPqe1tZWHHnoIYEUsFjt/uDGIlCsleERqXOJ0qjbSKmZOwBIqWbMtefDYtqRHgM3DiG8C1tgmV5+dfP0M68mT4h2B9zcN9nVC587Gdn2dlH7tIuACcseZJYa/Crz/9mBjEBERkeoQj8dPA54EmDdvHosWLaKxsbFgr7979242bNhAa2vPDKQb2ErfBbbkItt0sG1OS5cuzavPTj46Ozu59dZb8XaK6cRYLHZoqK8Vj8fHAd8B/gJg/vz5LFmyZMBkWGdnJ7fddhvd3d0As2Kx2PahxiBSzpTgERFC536AnQAF2JmSVzG8hEqqNmyL0gZsBpGPEViy5BwKs5f0j9jx6WmWA3dgh22Ny/M2Pv1FRmOlPCcOEMMjZGzT+m3g/Wvz/TOIiIhI9YnH408DCy+55BKmT59elPc4cOAAzzzzDJs3bybX57+xY8eydOlSZs6cOag+O/m4++672bt3L1jlzXeH0gA60RfoV8BpDQ0NnHPOOcyePTvv5z/wwAPs2LED4H2xWOy7g31/kUqgHjwiAnALKQmeQi9pzMA2PO/FkizrgK5+Hn86cDH599nJx6zsw+cnbkPWjCXD8jlHfh4ZCZ7LQuemBN7vGU4MIiIiUtFuAT7W2tpatATPhAkTOOecc1iyZEmy2TBAn2TP6NGjh9RnJx9TpkxJJni+AXwpHo8vA1ZiXQDyXWQbA9af6IILLsir8XSqmTNnJhM8bwSU4JGqpAoekQqVaEjcOdCJVHm+1mjswKmxybHrsCPQi+EothEc+lb0eKxyp5CJnaRDQKH3Qp0JvJz8M+Uem03s7zv80cD7rxQyLhEREakc8Xj8PGD56NGjufLKKwtePVMOjhw5wubNm9m2bRv79+8f+Ak5zJ49m3POOYeGhsHXKRw9epQ77rgD730XcFIsFstyyKlIZVMFTw0JnZsFnBR4/0CpY5GC+AAQhM7dgpWrPhjYL6xBC7xvD537LXBtcuwZipfgGZW4RWkccAP253qGjF44gzISO5lr0SCfl2xanVbF88HQuX8NvO/I9hwREREpL/F4fGwsFjtcwJd8FNja3t4+e+/evUyZkvUg0Yo2ZswYTjvtNE477TQOHTpEa2sr7e3tNDQ00NDQQH19fc/3uW719fXDSn6NGjWKlpYWtm7dWo+dEvaxgv0BRcqEKniqVOjcJOyU5vOwLSjnYQcRHQdmBN7vK2F4MkSJo7lvB+5x8B4Pp6Zc3omdBn4LcHfg/aD2NofOvQ34afL+ZOC9FK4PT7k5CDwL7MYSNiOx6qHk1xFZxpJfh3amhWkH/hPo7Dv854H3Px/Gy4qISJlKnPjTrKaulS8ej0+vr69/S1dX178B92Jzrl/HYrFhn4gZj8e/Dnxw4cKFnH766cN9Oclhz5493HPPPWCdA+YUOFEnUnKq4KkCie01Z2JJnGRC55QcDx8JXA38IJropMAuBV4DvCZLanY68L7E7UDo3O1YZc+dgff5nFZwO9ZseATYb73dQFMhoi5D44GzSvC+o4HFwOq+w38HKMEjIlIF4vH4DHoX15K3ifF4fEYsFlPPtcr2zq6urq8mvr8scfv3eDy+HJtz3RKLxTYM8bVvAT64bds2Fi9eXJXbtMrBlClTmDJlCnv27JkMvBP4VqljEikkVfBUmNC5emxnSHLicD6whMEl634feH95EcKTIgud+ynwtkE+7SjwHPAi1mcn21cHvAV4B5YEBOzY7wuHHbWkexH4YebwhYH3D0cdi4iIDF08Hp+AHfiYOi+bk+Ph74/FYt+JKjYprEQl1lPAwnPPPReA1tZW2tra6Orqs0P+CSxZ8ytgVSwWy+vDVjwen4ztIB/16le/etANhCV/W7du5ZFHHgFYDyyKxWLdJQ5JpGBUwVPmEgmdl2KHECVPjR4zzJd9ZejcjMD7tuHGJ9EJnZsCXJM6dibWrHgjtvcuh1EMvl0MYL1qlOApvCbsSPVNfYc/DCjBIyJSpuLxeCNwBn23v59K2m7mhoYGJk+e3HM7fPgw69atA3groARP5boIWDhq1CjmzJlDXV0dc+fOpbOzk7a2NlpbW9m+fTsdHR2LsWLdfwZ2xuPxVvpfZBuHLd69iUSLwtbWViV4imjWrFmMHj2a9vb2hVhl/B2ljkmkUJTgKUOhcw3YVpxrsO1UzYN9DQdMTTyxGXgMSKkJrsOa6X5z2MFKlN4DNCbvTAReif2/7gJewBIyzwJHCvSGnpQ9W1JQ55CR4HlL6NwXAu/XlSIeERHpKx6PjwBeQe8i2xmk/Uqsq6tj4sSJTJ48mSlTpjB58mTGjx/fZ3vN8ePHefLJJ+nu7n55PB6fqV48FesGgLlz51JX19uNr6GhgZaWFlpaWuju7mbXrl1s27aN1tZWjh07Nh3bQp+XpqYm5s6dS0tLS+Gjlx51dXWcfPLJycTrR1CCR6qItmiVidC5kdjn9WuwDP7UwTx/Ir3JnGZgBin7bIBHyDi554HA+4uHHLBEKnTuCuA2UiaWF2OzzXTdWH1v8rSoA4N8r7HYcuRiBjEjkUHzWCOstGYMx4GvAV/Is2+SiIgUUKJK59XYfOyN2JkDPcaPH9+TyJk8eTITJ06kvr5+wNd98MEH2b59O8BHYrHYvxUhdCmixDHmfwJGXXHFFYwbN27A53jvOXLkCMePH+fYsWM9X1O/P378OF1dXcyaNYs5c+YwduzYov9ZxBw/fpw77rgjub3uL4HvaKuWVAMleEoo0Rz5cqya5g1YnmZAo7Hjq1MTOgPt2doHfDdzeG7g/Qt5BywlETr3UuD3pPxvrgPej9X09scDh7GKnvbELfX71LEmLKkzl+GdEiX5W439j81iGxAA/xvoH2kRkaKKx+NjgCuxpM7rsT78gCV0WlpamD59OpMmTWLEiKHVtG7ZsoXly5cDPByLxbT7uYLE4/HFWHJnyty5c0n235HKt3r1ap599tnk3RXA38ZiseUlDElk2JTgiVjo3Djgtdgk4nVYwUS/RgLzgPlYYmcCQzu6+qdAWk3wPwXefzXrg6UshM6dgR3D2Sf591qG2FRHykoH8Eugn7NV7wU+FHi/JqKQRERqQqI58uuw+diVpCyiTJw4kdmzZzNr1qyC9UHp7Oxk2bJlyWqBebFY7PmCvLAUVTwePwm4H5jV3NzMhRde2Gd7llQ27z1bt25lzZo1HD16NDn8feDjsVhsZwlDExkyJXgiEDo3EavQuQZr5DVqoOeMAk7GzjqfS2GaJa0E7uk79Hjg/dkFeGkpgtC5Bdikos9OqcsA/U+rHh54Gsvk5NiT1Y31y4oF3u+NLDARkSqTOKXoKmw+djkpfe2mTJnCrFmzaGlpyWv7zVAsX76cLVu2AHwiFot9qShvIgUTj8dnYh0OTm5qauLiiy/OazueVJ7Ozk6eeuopnnnmGRKfjfcDMeA/YrFYR2mjExkcJXiKJHRuKrZ3+xpsL/eANb2jgQVYUmcOUOhfIYeA/8Q+UKY4NfB+fYHfSoYpdO5s7HjNuanjOra8eh3HjtB6FMvoZPEi8HHg+4H32iMuIpKHeDw+DetteC22RtKzZtbU1NST1BkzZrgHlA6stbWVhx56CGBNLBY7o+hvKEMWj8enAH8ElkyaNImXvexlQ96eJ5Xj4MGDrF69mra2noOGnwA+GIvF7unnaSJlRQmeAgqda6Z3EvFy8sjRjMOSOguA2RS/98nPsdOWUnwm8D5e5LeVPIR25MZrsN4rl6VfPwf7SzWU7XlSOfYAd5NxwlaqFcDfBt5rj7iISBbxeHwW8GZske1lJKZXzjmampp6TjwaNWrAguqC6u7uZtmyZXR0dAAsjsViT0YagOQlHo+/BPgGcO748eO59NJLaWxsHOhpUiW892zfvp01a9Zw+PDh5PDPgSAWi20pYWgieVGCZ5hC5+bQO4m4mDw+f0/AEjoLsZ46UX5gXwP8ru/Q08AiNXItndC5RuBtwD9gfY4zLMYyP0ru1AYPbMS2VO7P/bDvAx8PvNcecRGpefF4fC42F7sGeGly3DnH9OnTmT17NjNnziz5B/WVK1eyadMmgM/HYrFPlTQY6RGPx+uxdgoBVjDNmDFjuPTSSyOp7pLy09XVxYYNG1i/fn2yd9YR4IvAjbFY7Gj/zxYpHSV4hiB0bh69k4iX5POcydjWqwXYEeal+qDeDnyLjC0gVwTe/y7b44cqdK4OS3i9Hbga24GyDTvBO9fXfbWUaAqdmwJ8APgQdhhaVguwGYda+tWeTqxc55HE91kk94h/M/A+x0OKK3TuROAVWNXZ6Vhe6hbgwcD7rlLEJCK1IR6PL8DmYtdiha4A1NXVMWPGDGbPnk1zczMjR44sWYzpdu7cyX333QeWx18Qi8VqZt5TjuLx+GjgXcDfY1MuRowYwbx585g/f37kVV5Sfo4cOcKaNWvYtm1bcmgj8JFYLLashGGJ5KQET55C506lN6lzVj7PaaK3p04T5VN98Svgub5De4DzA+83Dve1Q+dOB67DKlJOGOTT2+mb8HkB6z37x8D7I8ONrVyEzp0EfAS4nn5OUZsNnIedoFYuf3ekNPZjPwgbcj+kFWvhszJ5C7x/sdBxhM6Nofev5mVYYuekHA/fBdyKHfWe4zR4EZH8xeNxhx0ieS02H1uSvFZfX8/MmTNpaWmhubmZhoZCHE9ReN577rjjjuSJPZ+OxWKfK/R7xOPxBuCV2HzsSuzXSK7FtVagNRaLtRc6jnKW6M30N4lbE1jFzvz58znxxBPVb0cy7Ny5k1WrVnHw4MHk0B1YoueZUsST+PfwVHrnYxOB3wC/jsVi/RzQKtVOCZ5+hM4twSYR15LnqdQz6E3qTCleaMOyCbg5c/gp4ILA+wODfb3QudnAW7GJRDGaBh7DPt/eCfwWWF+JlT6hc+djpb/XkKMgx2F/d87Ftu+JpNoM3IVlZPPwAr0Jn8ewpE/O7Vyhc6Ow5M2cLF+T308dYug3BN5/b4jPFZEalvgQcxa9i2wLk9caGhqYNWsWs2bNorm5uWJOONq+fTsPPvhg8u61sVjsl8N9zcR/p3OxudhfYFPSwdhL3+TPk9i864lqqjKKx+MLgb/DqnZGAUyaNIlTTjmFlpYWHYEu/eru7mbjxo08+eSTdHZ2AnRhHTBWptzWxGKxY4V+73g8PhFbPH8JvUmdXDsAVmDV1F8rRixS3pTgSZPYNnMd8F7gzHyeM5Pe7VeTihdaQT2Mnb+d5i7g7YH3OwZ6fqIh8CXAx4i+PcxmbNJxJ3BX4P3BAR5fMomtaq/HEjuX5HrcCGApdvz5xGhCkwrVBTwOPIjtexykrdjkYx22czQ1gdNUqBiz6AKuCry/o4jvISJVIpGsOJ/eSp2eSsGRI0f2JHWmT59eMUmddBs2bGDt2rVgfT0uisViq4byOvF4/GRs3nodNh0FYNy4cZxwwgm0tLTgnKO9vZ2jR4/S3t7e5/vk134+D2yjd4Htrlgstm8ocZZS4u/TRdhc7CoSc9bm5mZOOeUUmpqasGmtSH6OHj3KunXr2Lx5c7bLndg8Kz3pk7NvTzweH0/2hbXUr+PTn9fY2Mi0adOYPn06dXV1tLa20tbWluwZBPC/wNtisZhOX60hSvDQ8yH8MmzLzNXAgB34ZtOb1Mn4aasAHrgd67Cc5gDwaXL09Egkdl6HHdf80vTrqRqA+cBp2LL/oQFuQ2wg0onlqv4p8H7F0F6i8ELnRgPvwPZ0L8z1uLFYUucMEstIInk6DPwJK70r5W/teizJfQKWnNyEbQHNslx0BLg08P7R6KITkUoSj8dnAO/E5mM9vzsbGxt7jjOfNm1aVVRZeO9ZuXJl8gPiduADsVjstnyem0hYvBJbZHtlcryxsZE5c+YwZ84cJk+enHfSwnvPsWPH+iSA9uzZw44dOzh2rM+/5l3AQ1iy505gVTl/cEw0Tr4aS+y8BKw/0wknnMCCBQuYMGFCSeOTytfR0cH+/fvZu3cv+/btY+/evalbuFJ1YkeurwS2AC30TeYM+Jexvr6e0aNHM2HChJ6kzvjx4zN+zjs7O9mxYwcrV65MVhndGIvFgmH9QaWi1HSCJ3EC1nsStxP7e6zDPsCcgiUtcjZOqSAdWFo3R7nOGqx/xk6gLXGbDfwTVmySlQPmYkmdBUC+bQ099oHwEPbB9SA223keyzjl4Sjw3sD7/87zLQsuscXlZdh+9+uAabke24TVUZ+GfUAWGapOrNlNW8rtRYqT9HHAOCyR04L9mzgLq0BL1YUlnu7MfImdwIWB989lXhKRWpToF/MaLKnzemx9iMbGRmbPns3s2bOZOnVqVVZYdHV1cf/99/Piiz0t05YBH47FYln/jYzH43XAG7FFtvPAPvS1tLQwZ86cnlX8QvHes3//fnbs2EFbWxu7d+9Or/RpA74GfLWctnHF4/Fx2Nz+70hUf40cOZJ58+Zx8sknq3GyFFVnZ2dPsmeApE+PZPImeRszZkzG9yNGjBjUv4NtbW088MADyZ/Zj8RisX8b3p9MKkXNJXgSR1JfhU0kLqefrUXJZMVCLKkzOooAI3YI+CX2AXE4ZmBNik6lsMkvj20Kfz5x28qAlT5fBD4VeB/JilKiYfKVidtlQL9nac7FEjsnosbJUjydWJInNemzi/6TPsnkzfiUr8nbhMTYWAZ3mtvj2L7PNM8Crwq8z1rXLCK1IR6Pz8e2w7+bRNs55xzNzc2ceOKJNDc3V0WlzkC6u7t57rnneOKJJ5Kr7ceAn2Db0VP/GV+IVeycBpYAmz9/PvPmzYvslLCOjg527tzZk/Bpb+/py/wz4IZSNmpObHG5DHgt8BZsGzJjx45lwYIFzJ07t2wbb0v1SyZ99u3bx7Fjx/okc0aPHs3IkSOLksR+4YUXWLFiBdhHqj+LxWJZ2rBKtamZBE/idKfrsWO7++01MRE7lmExlbn9arC6sQ9iDzD4nh4nYTWvswsdVA4dWJLneWwrSI5ms3cB/4L15ynoX/C0Kp0r6Wf7VVIdlvg6F5heyGBEBqGL3qTPXixZk5rEGWzyJl9/ApZnDh8BvgDcGHiv5n8iVSx0rh77sH2pb2yc2nnaaed2z5jxMkaN6t2CNWoUTU1NTJ06tWZPL+ro6GDr1q3s2b2738eNbGxkxowZNDU1lTwBtm/vXp7ftInuri44duz5EStWfL1u9+59wMOB91m6AAxP6FwT8HJgnAe6W1pmdc+evbR74sSljB59CilF0WPHjaN5xgwmTZ5c6DBEKsqOHTvYtnUrQGfdpk3fa1i9+mFnH6cK/jlJykNVJ3hC5yZgJwlcjzXry6kB21K0BNsIWYvVFYewD2NP5vHYhdh/0MEe0VBo64DfYx9es3ga+Cbw46GcDpYUOjeP3oTOKxigSidpJNZb52xqI1Eoko3HzhF9KvvlZ4APBd5n2c0lIpUsdO5E4D0e3s2YMZuBM9yRI2p6Ujv+AHwDWBZ4n2Oa1r9Ej8xzsKqcKz2c72pzii5SDHdic7CSHPMuxVN1CZ5EE+CLsaTOWxjgw/gM4HSs3lU7cs12rC/PYaA98fVI4jYH2/RdTkfAtwK/xuLL4RDwI+Cn2B9vD3AwV9Z6KFU6SeOwqqaTsG1Y0RRNi5S3Luyszk25H3IL8HfatiVS2RK/P98EXO/hlT0fxp3DOYfvLtt+vFIkHjY7+BbWmHk3sCfwPudWrkSVzuVYQucK108/QxEZHg/HHYTAFwPvD5c6HimMqknwhM41A+/C9nOf0t9jR2EJnSVoy0y1OIB1JmzN/yld2E6VPYlb8vtJDKJKpw5rNptM6jShpSWRbLqBtcB9WEf0LNqxHlph4H3Oo0RFpPyEzp2BJXWuc+W1BiTl6Ri986/kHGwvcKqqdESi52GLs6bkv9K2rcpX0Qme0LkGrGzzeuzo7n4PJDoBS+osIHFEg1QVjyV4Hgc2ULyjo1OrdOYCjUV6H5Fq1A7cD6zO/ZCNwEeBuwPv90YTlYgMVujcJOCt2BzsnP4e24A1f18cQVxSWnsZ1GLboE1L3JQBEhmaPdh2hhx+D3wGWKkeiZWrIhM8oXMLsAnFu4Dm/h47HptQnI6VZkhtOISd874a22I2HKrSESm8HVg39H4mGWDJnkeBlYmvjwXe7y92bCKSXWIb/KVYtc61boDd7c3YwtqpaDGklrRhi21PM+DJpwMaiW15T87Bxg3z9URqncd+Nv9I7s9IHjqdFV6nzsHWBt4P9jweKYGKSfCEzo0FrsUSO5f099g67FjzJViFRfUfsim5dGFdXNdh/4gdxSoIBppwqEpHpPg89rP5J+znMk/P0HfC8fhwmqiLyMBC52Zhx5m/Fzi5v8eOAhZhczA1T6lt7dgnxI2J748mbgNVWE+jdw42iwHK80VkSI4BDwGPkd+uBw8dztbPU+dgTyjpU37KOsGTWCk6D0vqvJUBDiOaik0oFpFnAxWpWZ3YP2ztia+piZ8WVKUjEqV24AGs4m4Iv5E8tiszPelzqHARitSe0LkRwOuxap0r3QDrZSdi1dLz0TZ4yc0DHfQme1JvddjfI1XpiETnRayiessQnpto0ryavnOwJwPvOwoYogxSWSZ4Eh30346tFC3p77EjsdLfJVgpsD6Ui4hUpmRZ/w7sqJVh/HZKViCvxPo6/0L9fETyEzp3KpbUeacb4NEFqQsAACAASURBVCyKCVhSZzEwMYrgRESk4DywPnFrww6vGcZrHXOwCpuD/Q64QwmfaJVNgid0rh54FVat80YGOGG6BUvqnDLQA0VEpOJ0ALuwZE9b4usehpz0OQ78BvgR8H+B98NtCyFSVULnxgF/hs3BXtrfY+vpuw1eC2siItXlCDb3St52AAeH+Foedjn4GTYHW6VTuoqv5Ame0LkTgfdge7tP6O+xY7BVoiXoDE4RkVpznOxJn0FqA34KfD/w/okChidSURLb4C/AqnX+3A2wM2YaNv86DRgdQXwiIlI+DpOZ9BnCXvi1WKLnJ4H3OwsYnqQoSYIndG4U8CZspeiV9LMA5IB52KTiJNRoTUREeh0HdtKb9Gkj76RPF/D+wPvvFys2kXIUOjcdeAe2DX5Rf48diSV0lgAzULWOiIj0Okzf+dcO8ju92MMOB68NvH+8mPHVqkgTPKFzZ2BJnesYoAhnMr37utVsTURE8nUMS/psAZ4E9vX/8E8BX1DJsFSzxDb4K7BqnavcAH2QZ9O7DX5EBPGJiEh1OIQlep7HmiEey/E4DwcdvDnw/g+RBVcjip7gCZ2bhJ2AdT1wTn+PbcAmE0uwyYVWikREZDg8sA14AmsemOMsz28BHwy874osMJEIhM7NA97r4d3O2hfmNBZbWDsdW2QTEREZjk7gWWwOtonMPoqJo9ffHXj/s6hjq2ZFSfAk9nVfiiV1rgVG9ff4ZiypcyrQWPBoRERErHHzU9hxoFkyObcA1wXet0cblUhhhc6NBt6MzcFe0d9j6+i7Db7fc9BFRESG6DDwEHa8VhZB4P2NUcZTzQqa4Amdm4U1S34vcHJ/jx2FbfxegjXuExERicJWLJuTpWx4WeD9G6KOR6QQQufOxrZgvc3BpP4eO4XebfBjowhOREQEWAHcm/3ShwPvvx5pMFWq3z3Y+QidGwG8HlspupIBFoBOxCYV8wvx5iIiIoM0G9s3/Esyjv3sKEE4IkMWOjcZ62t4PXAm5N7e3oBVSi8BZvXzOBERkWI5D+uv+1ugu+8lzcEKZMg5ltC5U7EJxTuB6f09dgK9K0UTh/qGIiIiBdKENSR5uu/wg6WIRWSwQucuA27w8GY3wO72mfRugx8ZRXAiIiL9OBW4G0jbE685WIEMp4jmqf4u1mNVOkuAuWilSEREyktr5pAmF1Ip7oLcc6vR9G6Db4oqIhERkTzspW9yx8MhB+tKFU+1KfguqWnYhOI0bIIhIiJSbg4CB/oOHQceK0UsIoXgsG3wS7AmiPUljUZERCS79AU2Bw/rJNPCKUiCZySW0FkCzEDVOiIiUt5eyBxaGXh/NPpIRIZnIr3Hm48vcSwiIiIDyTIHeyD6KKrXsBM8E7Fjs0YMOxQREZHiOwL8KXNY27Ok4iwGXoMW1kREpDK8ADyZOaw5WAH1e+JVPg4AhTtoXUREpHg88HvgcN/hbuDnJQhHZFj2o+SOiIhUhqPY6VmpvOV8VMFTQMNO8HhgewECERERKbYngGcyh78YeL888mBEhmk70FnqIERERPLwB6wHYpIH7+CdgfeHcz1HBm/YCR6AbYV4ERERkSLaR+Loob4eBT4bdSwihdAFtJU6CBERkQE8BTydNuYgDLy/txTxVLOCJHieR9u0RESkfHVjZcEdfYfbgbcH3ndkeYpIRdhU6gBERET6cQCr3kmzGvhU1LHUgoIkeLZjWTkREZFytJys1ab/EHi/PvJgRApoBVadJiIiUm48tsB2rO/YMeC6wPtj2Z8lwzGcBE+fcqp7sJNJREREyskOsh7PcAfw7ahjESkEn5LT6cRWRlVJLSIi5eZRYEvamIOPBt4/UYp4asFwEjzvJyUZ1w78cbjRiIiIFFAHlsnp7jv8InB94L0+E0tFcvCPqfc3oUpqEREpLzuB+zOH/wD8e9Sx1JIhJ3gC7zcAn0sdexLYPNyIRERECuReYE/m8A2B9zsiD0akcL5Hlkrq9tLEIiIi0kcncDt2GECSh73AuwPvu7M/SwphuD14vgqsSx34HRlNLEVERCL3HLAqc/i7gfe/iTwYkQJKVJ+936uSWkREytB9wO60MQcfCLzXAdxFNqwET+D9ceAGUrZ+7ydrrwMREZHIHAHuzBzeCPxd1LGIFEPg/QaXVkn9BKqkFhGR0toMrMwc/kng/S8iD6YGDfsUrcD7R4BvpI49CrQN94VFRESGwGPVpGmN/7uwI9EPlSAkkWJRJbWIiJSNduzUrFTecj4fLEE4Nakgx6QDnwS2Ju8kj0M7XqAXFxERydcG4NnM4S8E3j8ceTAiRZSopH6fT6uk/mPJIhIRkVp2H5C6kubBO3hH4P3+UsVUawqS4Am8Pwj8derYi8D/oWM7RUQkWk9nDi0HPh95ICIRCLx/2ME3U8dWA2tLFI+IiNSmLmB92piDLwfe31eKeGpVoSp4CLy/DfhZ6th6YEWh3kBERGQAHcDzmcMfCrzXrhWpZp/waX/1/wBsL1EwIiJSe7aS0vkf8FbzEStRODWrYAmehA+Qtmh0H7CpwG8iIiKSzWbsaM4U27AKHpGqFXh/0MHVPuWk9C7gN8Dh0oUlIiI1JH17vIPbEluJJUIFTfAkmldejZ1xD9gWrWXAvkK+kYiISBZZeu+0JI6UFqlqgferHbw3dewQcCuW7BERESkWT9Y52NbMISm2QlfwEHi/EXgr0J0cO4qVCouIiBTTlixjoXPvijwQkRIIvP8f7GStHtuAx0sTjoiI1Ij9wMG0MQ9B6NziUsRTywqe4ElYT0qCB2Bhkd5IREQk6aLsw/8cOlcfbSQiJfNk6p06YH6JAhERkdowCZiRNuZgNDoePXLFSvB8CmhI3pkELCrSG4mIiIBtR1mZOdwJXBd4r10qUvVC50Z4+HTq2OnYPExERKRYtpLZksXDc8BHSxBOTWsY+CGDEzo3H+hTDn8hoKVTEREplueA35LSYbbXxwPv1WRZasW7HZyUvFMHvKSEwYiISHXrBh4BHsT68CR56HDwF4H3+0sTWe0qeIIHWznqyedMBk4rwpuIiIh0Yac1Ppr98m+Br0UYjkjJhM41eviUSxlbAkwsVUAiIlLVDgG3k73/oYOPBd6viDgkocAJntC5U4HrUsdeSvH2gYmISO3aB9wGtGW/fBfwjsD77uyXRarO9Q7mJO/UAxeUMBgREale/VROgzX7/9cIw5EUha7geRcp+ZxG1FxZREQK72ngd8DxzEtdWB+4Lyu5IzXm+tQ784DxJQpERESqU3+V0x52OXhn4P2dEYclKQqd4Pkp8LHknWPAJmySISIiMlzHgXuAtdkvvwC8NfD+wQhDEikXPwHOTt7ZhK2sji5VNCIiUlX2AcuAHdkv3+2scro1ypgkU0F3TwXerwN+lTr2EH0bLomIiAzFLuAmciZ3bgHOUnJHath/ediZvNMBPFbCYEREpHo8DfyYzOSOtz7L/wxcruROeShGe5zPp97Zji2pioiIDIUHVmMlonsyLx8D/ga4JvA+y2WR2hB4f8TBjaljj2E/ICIiIkPRAfwfVrmTvi3ewxYHlwbefyHwviv66CQb533h62tC55YBr0veHw+8GZhW8HcSEZFqdhTrtbMh++X12BGcqyIMSaRshc6N97DJwZTk2AnAG7G+iCIiIvnahSV2dme//Gvgei2ulZ9iHXD1udQ7B4H/IfsRaiIiIrncTM7kzg+Bc5XcEekVeH/QpZ1c8gI2BztUmpBERKQCHcIqp9OTO94KeT4IvFnJnfJUlAoegNC5m0g7Mr0eeC06WUtERAa2C/hR5vAh4K8C72+KOh6RShA6N9HDCgcLUscnANcAU0sTloiIVJBHgT+mjXnY4ODPtbhW3opVwQN2ZPp/pA50Abehpn8iIjKwpzOH1gNnK7kjklvg/X4HFwHLU8cPAP8NbCtJVCIiUkmeyhy62cE5Su6Uv6IleBKNlv4W+ET6tbuBe9HpWiIikp0n6+Tii4H3z0QejEiFCbzfBVyGtU/ocRT4BaAfIhERyWUP0JZy39u07COB99rtWwGKWcFD4L0PvP8S8G6gM/XaCuAOrKpHREQk1Xas4iDNpMgDEalQgfeHgauB76aOdwK3AlqCFRGRbLJUUHcCTZEHIkNS1ARPUuD9j4A3AIdTx58CfomO8BQRkb6mknIMUK8bQ+deGnkwIhUq8L4TeD/wmdRxD/wBuB9VU4uISF/zsN65SQ5GeLgldC7L1EzKTdGaLGcTOncOVrgzPXV8Gtb4b1xkkYiISLnbjZ3gcLzv8HZsD/j2EoQkUrFC597n4dsubXFvMXA5fSfzIiJS29YBd2YO3wm8PtGKRcpUJBU8SYH3K4ELSdv+vQu4CdgUZTAiIlLWpgJXZg7PBH4ROjcy6nhEKlng/XccvNFDe+r4E1hfnn2lCUtERMrQ6cAZmcOvAWJRxyKDE2mCByDw/jmynO5wCLgZ+D0Zq7UiIlKjFgDnZw5fBIRRxyJS6QLvlzl4hbcCuR5bgR9hfXm0ZUtERMA69c/MHP5U6NwbIg9G8hbpFq1UoXNjgf8BXp9+bQKWHjwh6qBERKTsdGP92jZnXrou8P5nUccjUulC507xcKeDk9KvzQWuwOZiIiJS2w4CPwGOpIx52O/g3MD7Z0sUlvSjZAkegNC5BuCLQAC49OtnA5cAIyKOq1g6gYZSByEiUoGOYFt5D2QOXxB4v7YEIYlUtNC5ZuDHwKvTr40EXoGV6GdMziqQx05t1RxMRGTwtgA/J6PCcy1wYeLERikjJU3wJIXOXQT8EJiffm0S1oOhJeKYhqoL28e+F9iT+Jr8/giwFHgZMKpUAYqIVKgdwH9j/86m2AJcHHj/QglCEqlooXMOeL+HGx2MTb9+ElbNUymHYByj77xrb8rNY4uGZ1KC/gQiIhVuJXBP5vDtwNWB9x1RxyO5lUWCByB0bgzwJeBD2a6fizVdKJdqHo992Gij7yRiPwPvXx+L7Wk8hepYGRMRicpa4P+yX/oTdrrDncCqoFx+uYlUgNC5k4DvAy9Pv9YIvBI4jfKZs3RhWzZ30zeZk88y8kzs1LBpRYtORKT6eGAZsD5z/IizqdlvgTsD77dEHpz0UTYJnqTQuUuBH5BlX/gUrJonS7OnyLQDTwJrSOtQOAQnY5Mm7XMXEcnf77B/g/uxA5ts3I19DmwFtgfeHyp2bCKVKnSuDvhrD192MCb9+nxsL1dGmU+E9mA/+0+QdhTYINUB52HHumrblohIfo4DP2XAz8BPYIttD2Lzr1ZgR+C9zlGKSNkleABC58YBXwH+Kv2aw05UifKXssf2AKwFNpCxPWBYRqCSYRGRwejEOvTvGPxTDwLbE7cDQEfidnyA7we63t/3hxOnR4pUhNC5+di2+YvSr43GFqZOjTCeTmzutQY77auQJmNJKx3qISKSnz1YT8QhZGtexOZfO7AcfTHmXOnfvxh43zaEP2ZFK8sET1Lo3KuA75Hld+9E4CXAYqC+SO9/GEtBrsVKfwdhKzYf2YBVsm0ANgJvBz5GltyUSoZFRPJ3GFgOPI9NNsrcywPv7y11ECL5Cp2rBz7s4QsuS9vAOdhC2xyKt21rFzb/ehI4mudzPHQ6eI7euVfythv4LPCmbM87HbgUS2CJiEj/tgOPY3Ow4VRTFpuH7Q4WB94P8qN8ZSvrBA9A6NwE4EbghmzXx2NltksoTH8ej9XzrwGexY7n7UcXtt/wEXonEc/01008dO504L+wuVEfdVh10gWoZFhEJF/7gU3YRGMztmxTZjYCSwPvjwz4SJEyEjp3KlbN85Js12dhc5aTKEyi5ziWmVmDfYDoT+KY3lvoLbBeD2zqr9ln6NzVHr7psuz2H4P1R1xI+fQaEhEpZx7rR7sJm4O1MnAv2hL4UeD9u0sdRJTKPsGTFDp3JfBdbD6RYQyW6DkDO95zsA4C67BZwoEBHov9Hf4u8MPA+9bBvldin/tfYU2lx6dfn4ydWjF7sC8sIlLjuoBt2GQj2XT1UOJrIbfXDsH3gA8G3pfzYpdIhtC5BuAfPHzW5ZhiTccSPQsYfHLEAzuxpM5T5FX2/wDwHeAXQ0mahs5NBP4F+Mts10/CKqozJmciItKvo9hC2xbs83RyDnaEkid+3gX8pFYO4KiYBA9A6NxkLCnyXnIU7IwCzgHOIvtR5B77C7cLm1Akv+4f+O07sJWi7wB3B94PUNwzsNC52cA3gDemX2vA9nM1DfdNREQEj008ksmeTizh053ytTvPsf4ek/y+nay/V3YB/wb8R62VC0vlC51bDHwdK3TJaipW6nMq2fsKdmFbKtPnYANlPT3sdvBj4LuB908OPvpMoXMXYxXVp6Vfmwa8jfI5uVVEpJJ1Y0meZLJnKHOtfOZeya8HyPp7ZQ3wZeDngfedRfhjlo2KSvAkhc6dAPwjtm0rWx6HkViSZz628XonvZOJY4N7uw1YUudHgfe7hhpzLqFzDrgaS/T0KRmegiV5hlKRJCIipXMc29eSoyL0EPBt4P8F3m+LLCiRAgiduxD4JPC6XI+ZhG05n0Lf+deLDLqS7m5sDnZL4P0gp28DC51rBD7q4ZPp1UmnA68p9BuKiEjRbQd+RvaqIQ+bHITAD6p163xFJniSQueagb8H/prCntx5FLgZm1TcF0U5V6Jk+CvA+1PHFwKvR/vBRUQqzWbgV/T7gbYDq0r4auD9+miiEimM0LmzsETPNYV8XQ9tDn4AfC/w/tlCvnYuiV5D3yetP+IVWI9HERGpLA8AD/Vz3cMuZ1Wp36y2quqKTvAkhc5NBT4MfAg7YGsouoHV2KTiplL8j05U8/wA2yfY4zLg7KiDERGRYTsErMR+ufTTW8RjW4C/HHi/PJLARAoksXXr4x7e6rLvzBqQh6MO/ohtmVrWX6PkYgmdm+BhhYNTkmP12FatGVEHIyIiw7YDO3F1Qz+P8XDI2e+efw283xpNZMVVFQmepEQVzF9jVT39ta85iO3DW4XNu1cBT5RDmVbo3Bgs4bg0OVYH/AU5ukuLiEjZO4r9slmJ7T/vxz3YHvHf1UozQKkOoXPzse1O73L9t6/ZTu/ca3Xi9kw59EQInTvdwyPOzu4AbNXwHeToByAiImVvL7ACeILcVdUeOhzcBHwl8P7pyIIrgqpK8CSFzo3Ftjr9I9ZLcxW9E4nHsWM0h90kuVhC5xYAjwITkmPjsQnGmFxPEhGRsteJTTCWM2Bz/1VYoufmcvjgK5KvlD6J78F2Kq5Kua0OvN9ZwvAGFDp3HTbJ73Ey8Ca0XV5EpJIdAh7Dfhnlqqr24B38GquqfiSy4AqoKhM8SaFzrlJXQEPnrsbaN/SYi210H1L9s4iIlI1urGR4OdaAth/PAZ8MvP+f4kclUjgVPgf7D+CvUscuwU4IExGRynYMS/I8hp2s2o97gY8E3q8qflSFU9UJnkoXOvdVIEgduwC4uDThiIhIgXmsxGE58EL/D7sq8H5ZNFGJ1LbE6Vr3AeclxxzwFuCEUgUlIiIFlayqXgHsy/EYDzscnB94vyW6yIZHCZ4yFjrXgB0Reknq+JuBeSWJSEREimU7luh5JvvlQ8BLA+/XRhiSSM0KnZvr4TFnp70Dtk3+ncC40oUlIiIF1o3NvZYDbdkfsgq4OPB+gIKf8qAET5kLnZuJ9Q3qOcRhFHADavgnIlKN9mCrSeuw0p0Um7FVpLLuYSJSLULnrvRwu0tpv3MicG3pQhIRkSLxWDX1w0CWcp1bgGvLuY9vktq5lLnA++3YIVo9f5mOYiu9IiJSfaYAVwCXZV6aC/wqsX1ERIos8P63Dj6fOvYCVtYvIiLVxWETrT8D5mdevhr4XLQRDY0SPBUg8P6PwI9Tx3aUJhQREYnIWcAZmcMXAf8VOqcDfUSiEQe2Ju90A7tKF4uIiBSZA14LTMu89InESYtlTQmeyvFQ6h0leEREqt9lZG3q+k7gl6FzS6KOR6TWBN53AX2OytUcTESkuo3ESnbGpI17+H7o3L+EzmXJ/5QHJXgqx4rUO5pciIhUv3rgKmBy5qWrgTWhczeHzi2NOCyRWqM5mIhIjZkAvAmbiyU5y/181MPzoXNfLsdEjxI8lWMd1n4HgMPAwdLFIiIiERmFZXNyNN65BlgdOvfL0LksO7pEpACU4BERqUGzsL6I6RyMBf7Jw6bQua+Ezk2POLSclOCpEIH3HdgRbT00wRARqQ1TsExOP7OHNwOrQud+FTp3ZkRhidSKlal3dgPHSxSIiIhEaxFwOZnbtQCcDf9joqLnq+WQ6FGCp7L0WUFqK1UUIiISuVnAO7Bqnhm5H3Y18Hjo3C1K9IgURuD9fmB96pjmYCIitWMp8D7gFfSb6AkSiZ4wdK6fqVpxKcFTWVQiLCJSwxxwMvB2Bkz0vAlL9Pw6dO6sSIITqW6ag4mI1LARwDn0JnrGZnlMItHzD4lEz42hc81RxghK8FSajMmFL1EgIiJSOoNI9LwReCx07jehcydGEJpItVKCR0REehI9N9Bvomc08PcenkskekZGFZ8SPJVlA3AgeecosKt0sYiISImlJnpe2v9DnwC2Fj8ikarVJ8GzFeguUSAiIlJ6qYmeE3I/rAN4IPA+stZtSvBUkMD7buDB1LG1JYpFRETKg8c68D+S/fJe4PWB958IvO+MMCyRarPK2yGmgH3zfAmDERGR0jsO/AF4IfvlVQ7OCbz/VZQxKcFTeX6YeudJLC0oIiK15ziwDLgL6Mq8vBw4K/D+9mijEqk+gfftDv43dWxNqYIREZGS2w38FCuRzuK/gJcG3j8bYUiAEjyV6NfY3ycAjgHPlC4WEREpkV3AT0g72qfX14FLAu83RxiSSLX7Tuqd54BDJQpERERK5yngJlI+lCd4OAK8I/D+A4H37dFHpgRPxQm8Pwb8KHVMK0giIrVlHbZqtDfz0kHgzwLvPxzlfm+RGvEI9uMH2PbIdbkfKyIiVaYT+D1wO1l30Tzl4LzA+5siDqsPJXgq03dT72wF9pQoEBERiU4HcGfilqWhzhrg3MD7X0QblUhtCLz3pFXxrEUnmoqI1IJ9wH8Dq7Nf/ilwfuD9kxGGlJUSPBUo8P4p4P7UMVXxiIhUv4fIWTHwPeCCwPsNUcYjUoNu8rZDHoD9gPZBiohUv1uBtrQxb+0Q/xLbllUWu3aV4KlcfVaQ1gEb0SqSiEg1W5R9+LOB9zeUaq+3SC0JvN/j4ObUsfvI7MMgIiLVZXGWMQeXBd7/Z6LCsywowVO5bsYWjgA4CtySGNxVqohERKSomoCTMofPDJ1zkQcjUrv6LLK1YUec3oV11xQRkeqzBGjMHD4z8kAGoARPhQq8P4IdoNLHZuDHwO+wkx3KJpUoIiIFcW7m0FXAjUryiETmT0Cf7ZAeeBzbK7mCrD2yRESkgo0ElqaNefjX0Lk3lSKeXFwZVRPJIIXOTQO+BLwXyDqxrwNGYH8hk19Tb7mujejne32CEBEpHY8l8rNUa34L+NvA++6IQxKpOaFzZwFfBy7O9Zh6Bj/H6u9xDWgOJiJSSgexEs7UiZaHLgfXBd7/b4nC6kMJnioQOncG8DXgsmK/Vz22RWAaMD1xm0bWcjURESmS7diW3GOZl57G2rJtSL0F3qtFiEiBJarmrvHwFZd192RhNdJ3/jUdmIrNzUREJBqPY1tyU3nodrCctPkX8Gzg/eEo41OCp0okJhmvB0LglKjffyK9yZ7kpGM8WmkSESmWncAvgDw7K++hd7LxTMr3TwfeHy1OhCK1IXSuEfigh085mBDle9dhSZ70hbfRUQYhIlJj1gL/l//Dt5KZ+EkmfwqejFGCp8qEzo3Ajmr7KPZ7fkSpYhlF74Qj+VUrTSIihfMiluQZxtLQs8CfB94/VqCQRGpWYuv8Z4C3exjnStjrcjyZC28T0cKbiEihPAXcwbB63i4D3hN4/2KBQgKU4Kl6iVWlcdjv+uTX8UMYS95GDSeeRuBlWIMqTTJERIZvL2nHKg5eB/CPwNfL6ZhPkUqWqKweTX7zqwHnZB4mOGvDM2RTgCuAluG8iIiI9HgGuJ2hN9b30OrgbYH39xYqJiV4ZFBC52YAZ2BHwiW/nsogV6lmA5djkw0RERkej23V2ovtxdqbcttH3hOPoqwkicjwhc7VYX1+zqTvHGzOYF/rDGyxTf0TRUSGzwMH6Dv3Ss7FDjBwhU+if8/ngM8H3g/7EEYleGTYQudGA6fTO9lITjzG9fe8euBC4Dy0bUtEpFg8dupD6oTjReCF7A9vpcArSSJSPKFzU7E5V88czMOigap9xgGvAuYXP0QRkZrViVVYpyZ/tgE5Tr74E3Ya19bhvKcSPFIUWVaazsRO+RqT/th5wJsjjU5ERNYDvyPrSVzdFHAlSUSildiefxq986+zgUuyPfayxEUREYlGF/AAduRWOg97nFVT3zrU11eCRyITOjcX+BZwZeq4A/6GYTb3ERGRQduH7R3fnv3ypwPvPxdlPCJSHKFzr/DwXy6taGcW8LYSxSQiUss2YU2aj2S/fGHg/cNDed2SdfeX2hN4vxl4HXAdVqEG2PaB50sVlIhIDZsE/AVwfvbLM6KMRUSKJ/D+HmdnXHwpdbyVYZ3CJyIiQ3Qi8C5gbvbL04f6ukrwSKQC733g/c+A/0wdf7ZE8YiI1Lp6rOHqyzIvjY86FhEpnsD79sD7TwCPpY4/V6J4RERq3VjgWizZk2bIczAleKRUfpN653lsP6KIiJTGpMyhCdFHISIR6DMH0yKbiEjpOLJOuIY8B1OCR0plObAzeec4sKV0sYiI1LyRmUOq4BGpTn2ad24GOkoUiIiIZJ2DKcEjlSXwvhu4LXVsY4liERGRwk4uRKSsrQZeSN7pTL0jIiKRK+QimxI8Ukp9VpCexRoui4hI9FTBI1IbAjtCN2MOJiIipdGYOaQKHqlIfwDak3cOArtKF4uISE3LkuCZGH0UIhKRPgmeEx22TAAAIABJREFUjWiRTUSkVAo5B1OCR0om8P4I8PvUMa0giYiUxmis0V+KmaFzQz6mU0TK2r0eDiTvHAG2lzAYEZFaNiZz6IyhvpYSPFJqGStIIiISvRFAlmxOltPTRaTSBd4fd/Db1DEtsomIlMastPselg71tZTgkVJbRkpVcBu2VUtERKI3J3Po5ZEHISJR0SKbiEgZGA00pdx3GUXV+VOCR0oq8L4NeDh1TBMMEZHSyJLguTT6KEQkIr/10JW8sxvYW8JgRERqWZY52JAowSPlQCc5iIiUgZbModOjj0JEohB4v9fBvaljWmQTESmNQiV4Ggr0OiLDcSvwpeSdLcBxsnYT79EFbMC2cx3BjuI6krh1AGOxs+WSt4mJr+OB+oKHLyJSHUZhfXh2ljoQEYnKrcBlyTvPAucO8ITDWCLoSJYb2FxrQpbbOIax50BEpMrNLtDrKMEj5eApbE4xHyx58zywMMeDD2CNe1r7ecE9/VxLJn9GYsmeeuwHIf379LHxwBQsWaTSNxGpVnNQgkekhtwK/L/knW3YotnoHA9+AbgdS/Lk8mKO8TpsLjWe7POuXHOwEcAkbA42FiWJRKQ6jQGmYttlh0MJHim5wHsfOncr8PfJsY1kT/BsxI58ODqM9ztM/xOTgdQBk7GJxpTE91MTX0cN43VFRMrBHGBlqYMQkUgE3j8fOrcWWAJ26sXzwKK0x3VjDRMfIuVkjEHqBvYnbkM1kt75V3IOlvyqDzUiUunmoASPVI8+CZ7nsEqe5HaqTuB+4NHo48rQjf3gZfvhG0PvpGMqNkHKtQomIlKOsvThEZHqdiuJBA/AM/RN8BzEFtdeiDiobI4DOxK3VA6rzk7OwZqxhUJVXItIJZkNrBrmayjBI+XiAWxn1RSwCp3fAK/Feu08RM7j0+9K3Ham3I5iPx9zU24nJL7OpIjVvck96FsT9x8CLgLOQJMMEakMo4FpwK5SByIiUbkV+GTyzjPY/OVMYAXwGLbQlsqDd/B9bJt9cv61C5vuzM1xm1ysP4Cntzro+cTYI8ArEm8sIlIJCtFo2Xk/1EJLkcIKnfsx8I7UsXpSzu/sqwubjHw18L57EO/RiCV/WrDPMY3YzqrGtO/Tx8ZhPYL+f3t3HidZWd1//HNmBnCYGRh2QVQQcJxhlzUCAsbgBoobKgRMjGFR+bndYH6Jxog/TNSrJlHzMwkRBEVRE2QxgoAOyiKyw6Agy7AN+2wMM8zWffLHqZ6+detWr1X31vJ9v1716q7n3q468hq7T53nec4zB9huzP+jarYikoydxvuDIiIlWkJ8mFvA8Ae6xF0tL0R6WGo2xeExi0mwDZrlYA7PGJyQuF85zveZRdRbtmX0vCv7/ZbAHIc5FvnYuOwKHE4bq0siIpPkxCrJW4idLDDx/EsFHukYqdnuwBWMvkPgMeC9ift17Y+qUWo2myj0vKr2GPp+V6IXYFO7AEegJENEOocTpxfezHBSkaUCj0jvS82OdzjXRsljgPlEcWeksy7aIjUzYAeG865sLvaykX52CrAfcDBRNRIR6QTriWWQt9K4cloFHukJqdkOxFLh/QourwH+FTgrcZ9s/6mWS82mATsTCccfAR8lDnyoMwV4NZFkqCmziFRlPXAPMVs00nYsFXhE+kNqdrjDf1ttu3yWw5MGZwL/nrg3WVxdndRsU+CVRA52DHBC0X2bAocCe6Ct8yJSnZXAHUS/nVVN7lGBR3pGajYD+C5wbG1oEDgX+Fzi3gk9/sakVqz6AvD+ouvTiSRjT5RkiEh5xpJU1CwCvpG4/2P7oxKRTpCa7ebwU4PdAByWGfwj8PXEfZRfGZ0jNTuIOP794KLr2xJb51vR70JEZKyeISbWfk/TNiRDbk7cD5jIe6jAIx0pNZsC/AOxq+kzifvvKw5pwlKzA4gk4zVF17chkowR1xaLiEzSeJIK4GvAjxL3dW0PTEQ6Smq2JTHRdgfwpcR9acUhTUgtl3yfwxetyfb/VwKvBWaXGpmI9BMntsDfwsinEToMGlxE5GDXJxMs1KjAI1KC2r7x9wJfIpo8N9iNaAKoJENEWsWJE2VuAR4e+dZB4CdEUnHdRJMKEZFOU1sZfobDGVawO34qsD9wELBx2cGJSM9aC/yO6K+zZIT7HFYYnE2slFw4wq1jogKPSIlqe8QT4K+JXVp1phLNhw5CTQBFZOLWAXczelIBrKCFSYWISKdKzV5GbDd7X9H1GcBhwO6AGo+JyEStAG4D7gRWj3zrQuCfgXMS9+da9f4q8IhUIDXbkUgymjYBHEoy1J9HRMZqBdFb5w7GlFT8C/DtViYVIiKdLjV7DfGhav+i69sBr2P0I11FRLKeJFZM30ssix7Br4kV05e0o2m9CjwiFUrNDiaSjAOLrm9LJBmFe7pERGrGkVRcSyQVF3fiSTgiImWo9ec50eEfDLYvuudVRH+ezUqNTES6ySBwP5GDLRrhPof1BhcCX0vcb2lnTCrwiFSslmScQKzo2aHontcRR6uLiAwZa1JBnIh+IfBPifvNbQ9MRKRLpGYzgf/r8Ekr2B0/jdjPtV3pkYlIJ1sD3EVshR9pGbTDEoNvAf+auI+SrrWGCjwiHaKWZHyK6NFT1wTwDcRx6iIiEFuwfgssH/m2JcC/Ad8sK6kQEelGqdnOxEEY78pfOwWYVXpEItKJ1hFLoe8imiiP4B7iFOXzE/dVbQ8sQwUekQ6Tmr2c2GmxYSbp/cRx6iIijwE/GPmWe4mk4ryykwoRkW6Wmp0InDf0fAZwKmq6LCLhOuCGkW+5ktgKf0XiPsqu+faYVsWbisiIBsgUd6YBW1UXi4h0EAd+2fzyVURScXlVSYWISJfbKPtkO1TcEZHwHHBTwbjDGoPvElvhF5QcVgMVeEQ6z2nZJ9uik7REJCwAnmocPodo2ndX2fGIiPSKWk/EU7Nj6r0jIkOuIZoaDnF4zuArBt9K3J+uKq48FXhEOkhqtgvRg2eDl1cUi4h0ljXEuZo5P0zcP1B6MCIivef9wAHZAeVgIgKxPf7e3JjBJxL3/6winpFoYYBIZ/kasPHQkxnkMg0R6Vu/AXINdVYDZ1QRi4hIL0nNNvc4zXSD3YAdK4pHRDpHk+3xtwLnlhzKmKjAI9IhUrM3Acdkxw4nU+0Rkb61lMgkcr6cuD9cejAiIr3nsxa74oHY4nBEdbGISAdpsj3+Y4n7QOnBjIEKPCIdIDXbmDj1ZoOXAHOrCUdEOsw1RPf1jEXAF6uIRUSkl6Rm8xxOz44dAGxeUTwi0jmabI+/MHEvGO4MKvCIdIYPAa/MDrwOndwg0m1WE6cseAtf82Hg/sbhMxL3lS18GxGRfpVapi/pZsCBFQYjIhPzHPBCi18zvz3eu2B7vJosi3SGd2Sf7I1ObhDpNk8AlwAriHN2typ4bM74CreDwC8ah28Avj+pYEVEhNRslsNR2d/LR5A7K11EOpoDdwFXE6udN6U4B9uU8eVgS4FbcmMGX0rcH5lszO2kAo9IxVKzGcDB2THNHLXHKuJD+EtRbyNpHQduJxrwDdbG1gFP1h5Z04AtGU42NmF4tc9g7fsdiC2aAHcAixvf8qOJeysXCYmI9KvDDKYOPdmCaK4srbeUWGHxMrRCXVpnHXAVcHdmbFXt8Wju3ukM519bEP/H99xjLjCzdv98hvM64vpjBl9q7f+C1lOBR6R6h5CZLJqN9n230jJie8v9RNMSJ4o7uxMrpbauLjTpAWuBK4Hfj/H+9cDTtUczmwJ/RiTA1zVePjdxv2k8MYqISFOvyz5R8aF1nJhUG8rBltTGNyfyrz2Iv3ciE7WUWDn9zBjvf4E47vyxEe55GHhn7esDuWvWJdvjVeARqV5DciGT8wxwD5FQFKx+YC1wW+3xUiLR2I3MFJ7IGCwBLqb439hkrCKKRjOJjd4ZzwN/0+K3ExHpZ8rBWsiBhcB9xIfjVQX3LAd+RUxgzAH2AbZHhTUZn/uAnxE5fSs9RKyevq3x0vXAD1r8dm2hAo9I9eqSi4XEH8Vdqomlq60FrqXwOOmmHq09ZgB7EsWeWa0PTXrMvcAVFCYWzwMfIP4pzss8dq89thzL699XPHxW4v7E+KMVEZG81GxLh32yhYXbiJW9W1UUUzdbSkxOjLU5yQDwu9pjG6LQMw/1P5KRDRKnWhUtZXa4x+A4YrHO7mTyMIe5Fru0RnVV8XDXbI+3LolTpCelZrOIXUQNJ9rtCRyFZjTG6gHiF/KK0W9dSdRzCh0NvKpVQUnPGSBmHvNN92ruBt6ZuN9bdDE1MyKPHUo45hA7BgeInOX1xPbvIguBeYn76ibXRURkHFKztwE/yY9PAd5AfDqU0Q0ANxPd/9ePcJ/H5UFr0gZxCnAKIyRo0vdWApfSdIvVD4EPJu6FHwVSsynAyxmedNuZWHQ2CLjDn1vzf37nJu5/PqngS6QCj0jFagnG14ndQnX2Jj7xqcjT3ErilKHCT9Rh6DP5xcRW3UXEqWWnAa/N3jidSC60tFGKrAAuI/4BFfgecMpk9manZvOIBWibFFz+y8T97Im+toiI1EvNNgYSh88YvCh7zYC3oAmf0TwB/JzmPVAcVljspLmY+DqVaDN3GvCK7L27AW9rW6TS7R4jijv5JMthvcEnga9PZoVNanYK8K0ml1/e6SdnZanAI9IBUrOZwGeBj5NrBbMfcWSnijz1ho5EvAZYU3zLQ8DfA5cm7kuKbkjNfgscMPR8f+K/tUjeI0Rxp6CfwDrgo8C3WrF0NzX7BPCVgksvTdxH6gsoIiITkJq9Avgm8Mbs+BTgrcCuVQTV4cawJf5q4MvA/MS9IU1LzbZ2eMIyc2rvAnZqdaDS9ZxYIfYrhk8dzVxbZHBc4n79ZN+ntsr6f8j9HgB+n7jPm+zrl0kFHpEOkprtQ/xy2T47fhBwWCURdabFxD7vJp92B4CvAp8baTVFarY7sCA79gHG2CBF+oYDvyUS2YK/lo8A707cf9uq96stIf4l9avL/pC4z2nVe4iISL3ah7uTHP4ze2z6VOBYYi+HhJG2xDssNvgEcP5Ikx6p2ceArw093xz4IJrMlHprgMtp2pfwF8D7EveRDiYdl9TsJQ4LLA41HvKNxP30Vr1HGRr6fohIdRL324E/Jrfa9UZib7PEf4vzaFrcuQU4IHEfyzGGJ2efvBQVd6TeaqI5w68pLO5cAezXyuIOQOI+CLwH+CeG8+dftPI9RESkXuLuift3LIo8G37lDxB7i7pmb0YbrSH2uV9E036H3zWYm7ifN0pxx8jlYHuh4o7Uewb4Ls0PnQCOamVxByBxX2TRAvUij9480IU5mFbwiHSg1GwvYD6wRXb8CGIbUb+6n4JuiGEV8Gli/+1IPf4ASM2mA4+TqdCrubJkPU0k9csbLznwOeD/Je4D7Yyh1oT9z4C7E/euSzBERLpRavYB4D+zYxsB7wZ2qCSiznAVcHvBuMNCg1MT95+P5XVSs8OIHTeAmitLo7uJlfr5hN5hmcGJiftl7Y4hNdsJ+DDwhcR9abvfr5VU4BHpUKnZ/sQe5s2GxnYi9ij3o7XAORTOGv0M+FDi/tBYXys1Own4ztBzNVeWrHuIJcEFlcIlwPGJ+xUlhyQiIiVKzT4MfCM7djiZpn195nHggtyYw4CNYUt8Xmp2PvCnQ8/VXFmGOPHBp6iQSLR8elfivrDEkLqStmiJdKjE/WZiN9IG21UUSye4lobizgBwAvCW8RR3auqWBu+OijsSniUqhgXFnZuAV6u4IyLS+xL3bwJ/yI71aw42QJyUleXwsI19S/wGqdmWHouhNti7FUFKT7iDpsWd/wAOUXFnbPSZRqSz1f3d69fk4kngtsbhNHHPTyiNqtZc+ZDs2F4TDUx6ynrgp0Qym/OvwCeKTgIREZHek5rNcNg12xemX3Owm4nJjyyDv0jcC1KzUZ1ksMnQk82Bl08mOOkZi4neFFkOqw1OS9zPLT+i7qUCj0iHSs2mAq/OjvVjcjFIzBzlNpMuBM6c4EuqubIUupZcd/NwWuL+rdKDERGRKu1jmZ0OW5CpSvSRZRQe8nF+4n71eF9LzZWlmQHiCOHs6mmH5w1eO8FCYl/TFi2RzjWHTM+56WSa8fSRW4iGtzmnJe6rxvtatebKJ2XHtDRYAB4mZilzLlBxR0SkL9WdadGPE2xONFbOfeheDHxygi95KDB36MkUYI8JRye95DrgqdyYwekq7kyMCjwiHSg124E4HXCD7ei/WY7lxC/9nO9Pog/Ku8mcnDUd2HWCLyS94wWi707OI8TpCSIi0kdSs8MdPpsde3FVwVToHuCh3JhBkrgXLHYdk7rVO7ugk7MEHgV+2zj8YzKHocj4qMAj0mFSs72BG4F9s+P9OHt0Mw3NbpcBH5/ES6q5stRx4ijO5xuHT0zcl1UQkoiIVCQ1O9HhSotdWRv0Ww7mFG7Nms8EP3SrubIUWU1szcpyWASckuio7wnTZxuRDpKavRH4ETAzO74tuWpPH1gDLGgc/uvEPb+Kc0zUXFmK3E3umJTwxcT9V6UHIyIilaj1h/ks8Nn8auk9gR3KD6lSDwNLMs8d1hucOokP3WquLA2upuGEXAzen7gvKbpfxkYFHpEOkZodB1wATM2OvwI4Gti4iqAqdDewrn7oSeCcSbykmitLnWVEcpFzK7ml+SIi0vP+P3BKfvAw4ED6b4v8rbnnBj9K3O+dyGupubIU+X3tkfPViTTwlnoq8Ih0gNRsc+Ab5Io7+wJH0n97KZ3CY9G/lbivncjrqbmy5A0Sy4JzRcQXgBMm+u9MRES6T2p2BLnizlTgTcCrqgioYsuABxuHvz6Jl1RzZamznNgen3Mn8Ddlx9KLVOAR6Qx/C2yTHTgS2K+aWCr3ELC0fmgd8G+TeMl3oebKknEL8Hjj8CcT93tKD0ZERCqRmk0Fvpodmw4cC7ykkoiqVzDBdgvwm0m8pJorywZDvQ/X1o+tMTg+cV9TUVg9RQUekYqlZq8APpod+yP6s7jzLNF35+7GSz9K3J+cxEvXzcypuXJ/W0Vh88jLAB2JLiLSX04k1+awH4s7g8TRkQso7Ev39Yn23hlqrpzdjqUV1P3tQQpPZzsjcS9I/2Ui9BlHpHoJmRY7M4EDqouldKuJozgXEE12mviXib6+mitL3nXUzxwRq4U/qBMbRET6zqezT+bQX8WdZQxPrOWb3QI4PGtw4STeQs2VZYMB4ii2nPlEmwppERV4RKo3K/tkHr3fUNmJExoWAPfTcBR63k2J+42TeDs1V5YNniE2eeecOdHT2UREpKvV5WD7VxVFidYSq3QWAI+Ncq/BvyfuqyfyPmquLHm3U9+CwcENPpa4D1YVUy9SgUekenX1jfuJUxt68Q/gaDNFOc8BPwD+eaLvp+bKkuXAL2tfM+5HM0ciIv2qLgd7BNi+okDayYFFRA52Lw0HDBR5CDgXNVeWFnkBuD43ZnB24n5HFfH0MhV4RKpXt1tkCXAe8Gri9IaNqoiohcYzU1RzNXEc+kWJ+6pJvr2aK8sGDxLJe84ndWqWiEjfqvv9/2uiAf++xFaibp9sW0FMqt1Nw+EVDRxeMPgxkYNd04JVFWquLBtcD2Q7KDusMPhMVfH0MhV4RKr3eeAoYKehgWeAK4BrgB2IwsTQY9Pc8+nAi6gmCRkgGtY2e6wkZozGMVP0ncT9oRaGqObKAsQU7fzG4auBS0sORUREOsepDpdaZj7tgdpjNrA1xXlXNh/biGpysLVErvUCjbnXKqK48zgNq1aL3AB82+CHiftzrYhNzZUl61lie1aWwVnaHt8epp6SItWrnaT1a6KeM25G88Sj2VjRyiAnquvNijX5RGJCm7KHvUBrZ4rq1JorL8iOnQhs18o3ka6wDriYhlMbBoF9E/eCljwiItIvUrNjHX5sMHUiPz+V5kWgfO61KTEpV/RGgzQWa0Z6jGHyrCmHJywWjJ+buN8ziZcqlJp9DPja0PPpwIfo/hVRMn5LgB9R35rBYaHBvIn2d5KRqcAj0iFSsznAF4FjiK3KbTWN4WTDGU4YSuhydgPwbVo4U1QkNfsq8PHs2LbAO4iTyqQ/rAH+m1hJlvPvifspjcMiItJvUrO3ECuq9x3t3lbYhOEV2OsYnkBrJ4d1FvMd5wA/T9xHOeNi4lKzBcTC6Q32AP6ECVbRpCs9QxR3CvotvDtx/3HZ8fQLFXhEOkxqthNwGvBBeufApydo40xRkdRsLvA/ZLa+QRyX8U5i2bX0tlXAfwEF63//AByauD9TckgiItKhaqc+/RHwkdr2ol7Z1X07UdS5IHF/tow3TM2Oczgve0Q6REJ2DLlB6UlPEMv01zReuhB4X6IiRNuowCPSoWonQO1H1CKGHts0+X5Wk5dpNye21j4FPJ17DI09CdzazpmiZlKzFwOXEf8dN9gEeBvwsrIDktI8T8waLW68dAdwVOL+dMkhiYhIl0jNtgPmMXLutbXDNvkiRllqK3Ky+VZRDvZQWRNreanZoQ4XW26ychtiNXVViau03yPARRRuI/w2cHLiPlBySH1FBR6RHpCavQjYioLkg+aJSbMDulbSmCAUJQ1PA4s7/Zd0ajYT+D5wdHZ8CvBGInuT3rKMKO4sb7x0A/CWxH20w0RERERGVVv1M4OR86267x22tObtaJYy9hxseaevgkjN5jj8zGDn7PgsosizTTVhSRs9AFxCHMSS88/AJ1rdc1MaqcAj0odqCclmDCccRi1pSNxXVhlbO6Rm04CvA6fmrx0KHIQa//WKVcRewOcbL10NHJu4F1wSEREpR2o2FdiCqG9sRUysPQU8m7ivHelnu1FtNdSlwAHZ8Y2BY9Fq6l7yMLE1vqCCcybw951ekOwVKvCISF+oFbXOAP4xf+0twNzSI5J2+AMxc5RzCfAendYgIiJSvtRsBnAB8Nbs+BTgZHT4Ra+4ktgHn/NXiXtaejB9rO0n9YiIdILE3RP3LwLHA3UzZEuqCUnaoKCh8s+Ad6m4IyIiUo3a6vB3AN/Mjg9SuJ1aulRBDnaGijvlU4FHRPpK4v594G+yYysqikVaryC5+EHiXtDnT0RERMpS69l4OnBFdlw5WG8YII5Fzzm/9EBEBR4R6Ut/yD5RctEbnMICzy2lByIiIiINaj1YlIP1oGdpaKz8eOL+ZCXB9DkVeESkHz2afaLkojesAF6oH1oFVHI8rIiIiBR6LPtEOVhv0ARb51CBR0T6UUOBR+3mu19BcnFbbUm4iIiIdAZNsvWgghzs5vKjEFCBR0T60xIyiz3WAWuqi0VaZGnj0LPlRyEiIiIjqCvwPFdVFNJSyxqHlINVRAUeEek7tT3gWiLcY3ZsHDoyNdu4/EhERESkCa3g6UEvaRw6qvwoBFTgEZH+pRmkHrMNDX/UNgMOqyIWERERKfS4Z3bGrwLWVxiMtMb2uecOf5KavaiSYPqcCjwi0q9U4OkxNwGDjcOrSw9ERERECiXu6wyeyI5pFU93c+A3jcPPARuVHYuowCMi/evB7JNHm90lXWExcGPj8H8k7teVHoyIiIiMZGH2iXKw7nYnsCg3ZnB64q7aXQVU4BGRfnVF9slCtES4m90A5I7Legr4VBWxiIiIyIjqcrD7qopCWuLaxqHLgP8qPRABVOARkf51E5klwuuAh6uLRSapYB/W3ybuBQdriYiISMV+kn3yCLC2okBkcpzMsbTDPlI70EQqoAKPiPSlxH0QuDg7dn9FscjkFay+0oSgiIhIZ1pAZqv8ALk9W9I1BhqH1ibumjOtkAo8ItLPLso+eYDCJr3SBQoKPK9Nzaz8SERERGQktdUddTmYJtm6U77A4zA1NTukkmAEUIFHRPrbfDIHaK0ilgmvAp6vXVgGLAGerX1Vn57OVDCD9HngvNRs09KDERERkdHUbdN6kMi/8jnYYuAZYDmahOtE+bzYYKrD/NTsdE20VcO0PU5E+llqdgHwvvH8zGbAbGCL3NfZwLRWByhjcgPQ5LisO4F3JO4PlBmPiIiINJeaTXV4wmCbsf7MVGBz6vOvoe9noZULVVhHVOqa7Mn6HnBy4r6qxJD6ngo8ItLXUrPjgAtb9XqzgC2BXYDdgU1a9cIyqvuAn1HYqHEZ8KeJ+09LDklERESaSM3OBv6iFa81VPzZCphH5GEq+JRjkJhku7H4sibaSqYCj4j0tdRsM2L178atfu1pwKuAvYEXA1qn2n5LiM7Zi4svfw44s9ZgW0RERCqUmh0NXNqO154J7Fl7bNaON5AGzSbaHJaZJtpKowKPiPS91Owy4Chipen6Eb7OBHZkArWabYG9iFmllleSpM5a4Arg3uLL/wOcmLgvKTEkERERyUnNXgQ8BbyI4txr6Pv1RCq19Xjfw4Cdicm2ndGqnnbTRFv1VOARERmHWjLyCmA3YNfc15cySvFnI2I26WBA3X/bx4FbgGtq3+csBA5J3J8oNyoRERGZqNRsNo2519DXUYs/s4ADiGLP1PaF2fdGmWi7GHh7oiJE26jAIyLSIrXiz67A24G/JAo+hTYmijyvRo2Z2+kR4DLiVI6c8xP3k8qOR0RERFqvVvzZAzjJ4XiDGc3u3QI4nOjTo+3z7THKRNuJift3Sw6pb6jAIyLSBqnZVOBNwCnAm2myKngz4DCiV4+SjPZYAVwCFCzXOSRxv77seERERKR9av0VjwdOJRbsFNoROILokyjtUTTRVjs9bU7ivqKisHqaCjwiIm2Wmr0M+GDtsX3RPS8mkowdywurr6wnzup8pn74VuDAxH2ggpBERESkjVIzAw4ETnF4r8H0ovvmEpNtasbcHkuBc4FcsvXlxP2MCsLpeSrwiIiUJDXbBPgI8GlgdtE9uxFJxpYlxtUvHgUubBz++8T9c6UHIyIiIqVJzbYHPufwF1awqnoqsB9wELBJ2cH1gWuB32SeO6w3eH3ifk1VMfUqFXhEREqWmm0B7Af2AAAJ/0lEQVQFfAb4ME1a8LyUmFF6JXG0hLTGpRQ2/fuzxP07pQcjIiIipUrN9gC+DLyx6Po0opniXGAn1Iy5VdYC5xDb5ofUjk8/JHH/XUVh9SQVeEREKpKa7Qb8A/DOZvdMJY71nEs0A1RD5slZAXwHWF0/vB54U+J+VQUhiYiISMlSs6OIQs9eze6ZTky0zQN2QL0SJ+s+4gitLIdHDA7WyaatowKPiEjFUrNDgK8QK4Ob2phINOYSK3wKuzbLqB4DfkTDXvAVwKGJ+50VhCQiIiIlqx2IcZLDWdakR+KQzYj8ay5jOI9dmroR+HXj8G3A4Wq63Boq8IiIdIBaI8D3AJ8ADhjt/hnEyVtzge3QrNJ43Uts18pZRMwiPVZ2PCIiIlKN1GwG8HGHky3m0Ea0DZF/vQo1Zh4vB64ECmbTLgeOSdzXlxxSz1GBR0Skw6RmrySO9zyB2Ao+ollEsjH02BrYAu0bH83NwPzG4buAwxL35SWHIyIiIhVKzaYAhwInOLzbIp0a0VbU51/bEHmZJt6aGwQuAhY2XjobODlRgWJSVOAREelQtVU9+xOFnvcSi3XGZCpxEtdQwrEbY8hS+owDvyDWBedcDbw5cV9bckgiIiLSAVKzjYlGzCc4vNXGcebFJkTuNVTwmUdss5dha4mTTZ9qvPTpxP2ssuPpJSrwiIh0gdRsGvA6YmXPO4gJojGbSVSJxvVDfWAQuAS4v/HSecTpWvojKSIi0sdSs1nA24liz+uLjlkfyc61H1bvxHorge8BzzVeOilxP7/seHqFCjwiIl0mNZsOHEMUe94MbDSWn9uWWAakWaR664AfAgXHN3w+cf+7suMRERGRzpSabUf0TDwBOHCsP7cP8Mdo61beYuD71J9u6rDe4I2J+9UVhdXVVOAREeliqdkmRJ+/vYA9M193KLp/F+BtaBYpbxVwAbCs8dJfJu5nlx2PiIiIdLbUbDawB5n8y2Eva7Jg+khgvxLj6xaPAj+m/nRTh+csTje9q6KwupYKPCIiPSg124pIOD5F7CHfYE9ir9eYlv30kaVEkeeF+uEB4OjE/fIKQhIREZEuUuuf+DJgb4dv5E/legNREdJKnnr3AJflxhweszjddFEVMXUrFXhERHpYajYT+BWwb3Z8M6LIswtKMrIeJ7Zr5c7ofB54beJe0I9ZREREpFFqtpfDtfkVPTsCryeaMMuw3xIJa86dxOmmBa16pIgKPCIiPS412wG4kcgp6uxM7AmfXXZQHew+4OLG4SeBfRL3ggMfRERERBqlZm9w+KnFAacbTAFeDbwG9UYc4sQxprc3XrqcON1UhYsxUBsGEZEel7g/DhxNQR/hhcA5wHVEs2GB6RSuarqP6AUoIiIiMiaJ+xUGJ3ucDL7BIHAz8G3gXqK40e8M2LT40o0q7oydCjwiIn0gcb+DaMb8Ver72DEA3ACcCzxSemSdZRWxBzyXRSwGjk/c1xf8iIiIiEhTifu3LVrvXJG/9jxwKdFkuN/3ID0CXN84/Evg82XH0s20RUtEpM+kZnsC3wQOy18z4Ahi2XC/9eZx4L+JVU05b07cf1Z2PCIiItI7ag2Y3+7wT/nmyxCrV95KwX76PrASOK/2dYjD0xbb4xtWoEtzKvCIiPShWpLxp0AKbJu/Pg/4E/rrpK27gYIqzpcS90+VHoyIiIj0pNRsBvAZh08aTMtem0IcgrE3/TXR9lPg95nnDm7whsT9yqpi6lYq8IiI9LHUbDax9PXD5HKJ7YC3ESdu9bp1xD74FfXD1wNHJO5qTyQiIiItlZrNJVZUH5m/tidxCMa0/IUe9CTw3cbhsxL3T5ceTA9QgUdEREjNjga+R66e0y/LhW8Efl0/tBaYm7g/WEU8IiIi0vtSsynAZ4G/y1/bnphom1l2UCVy4IfAo/XDdwL7qffhxKjAIyIiAKRmc4gTwudkx3t9ufAq4Gxyx1vAVxL3pIp4REREpL+kZm93OM9y9ZwZRJFnh2rCarsHif6HOW9I3H9eejA9QqdoiYgIAIn7vcBBxIEOGwwCVwE/B3ptKsWJlTu54s4y4AsVhCMiIiJ9KHG/yOBgh/uz4yuBC4G7qgmrrdYC8xuHr1RxZ3JU4BERkQ0S9+XAscCZ+Wt3EUnG82UH1UY3UZg0nZW4Lyk9GBEREelbifvdBgcCl2fHB4jz1a+qfd8LBoBLgGyy5THvdkY1EfUObdESEZFCqdnbiVMre3K58AJyGVR4iOi9s7rkcERERERIzaYCZwENp3juCBxD5GLdyolTS3/XeOm8xP39ZcfTa1TgERGRplKz3YGfALtmx6cSe7kOoPuOUh8Ebgd+SSQZGc8Bhyfut5cflYiIiMiw1Ow9DucYTM+OzySO3Xol3dcbcTWRf93deOlW4MjE/bmSQ+o5KvCIiMiIUrMtgAuAN+avbU40YN6l7KBqnEgWVhNJjhF7j/PfD319Gri69jVnLdHUb377oxYREREZXWq2t8NPDHbKX3s5kYNtVXpUYZA4qGIdxTlX/vt7gF8BL+Rex+EBg0MS96dKCr2nqcAjIiKjGmm5MESB50hgdgvf04nEYQXR92dFk+9b0PjZgeMS9x9P/qVEREREWic125pog/i6/LUpwP7AwcDGLXzPAaLBc7Pca+jrZCsJDk8bvCZxf2CSLyU1KvCIiMiYpWbvBr4JbJO/No3oDHhg7fvxWA48WnssZzhxKKGZ4Arg9MT9O+1/KxEREZHxS82mAWc6/JUVpFmzgCMY/7YtJ1Y1PwosYriIs3KyAY/tve+3mGC7rYS36xsq8IiIyLikZrOBzwMfouA0xrFs21pBJBOPMFzUqcD3gSRxf7yatxcREREZu9RsLvANClbzwOjbthx4luEc7DFim3uZHFYbfAH4sg61aD0VeEREZEJSs32I1TyvKbq+C7Azw8t3s0nF0taGshJ4pvb91MxjWu750OM+4COJ+y9aG4aIiIhIe6VmBhzn8FUrONR0CrAvsAXDOdh64AkiB8v3wJmkxcQ8XVG+VfS4HPg/ifvC1oYhQ1TgERGRCUvNpgAnAV+iYNtWCywnJphGeixP9MdMRERE+khqNgv4jMPHi7ZttcBTjJx/LUrcW1wvkslSgUdERCZttG1bY7AW+A1xeuYNwMNE4rCiZUGKiIiI9JjRtm2NxmGZwXwiB7uNWOjzeOK+tmVBSmlU4BERERERERER6XITmWUVEREREREREZEOogKPiIiIiIiIiEiXU4FHRERERERERKTLqcAjIiIiIiIiItLlVOAREREREREREelyKvCIiIiIiIiIiHQ5FXhERERERERERLqcCjwiIiIiIiIiIl1OBR4RERERERERkS6nAo+IiIiIiIiISJdTgUdEREREREREpMupwCMiIiIiIiIi0uVU4BERERERERER6XL/C4VcySpHzdsrAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 145,
       "width": 572
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(ncols=2, figsize=(10, 5))\n",
    "patch = PolygonPatch(\n",
    "    world.loc['Belgium', 'geometry'],\n",
    "    facecolor='red',\n",
    "    edgecolor='darkred',\n",
    "    linewidth=2\n",
    ")\n",
    "axes[0].add_patch(patch)\n",
    "set_limits(axes[0], *world.loc['Belgium', 'geometry'].bounds)\n",
    "axes[0].axis('off')\n",
    "\n",
    "\n",
    "patch = PolygonPatch(\n",
    "    world.loc['Belgium', 'geometry'],\n",
    "    facecolor='lightgrey',\n",
    "    edgecolor='grey',\n",
    "    zorder=1\n",
    ")\n",
    "axes[1].add_patch(patch)\n",
    "\n",
    "redpatch = PolygonPatch(\n",
    "    clipped_belgium,\n",
    "    facecolor='red',\n",
    "    edgecolor='darkred',\n",
    "    linewidth=2,\n",
    "    zorder=2\n",
    ")\n",
    "axes[1].add_patch(redpatch)\n",
    "\n",
    "set_limits(axes[1], *world.loc['Belgium', 'geometry'].bounds)\n",
    "axes[1].axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now let's try some countries we might be a little more familiar with (with a couple of geometry tweeks to make it look nicer):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "world = world.loc[world['TYPE'].isin(['Country', 'Sovereign country'])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dirty hack for assigning a multipolygon - see https://github.com/pandas-dev/pandas/issues/26333#issuecomment-561846798\n",
    "\n",
    "world.loc[['New Zealand'], 'geometry'] = gpd.GeoDataFrame(\n",
    "    geometry=[ops.unary_union([p for p in world.loc['New Zealand', 'geometry'] if p.bounds[0] > 0])]\n",
    ").geometry.values\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "world.loc[['Netherlands'], 'geometry'] = gpd.GeoDataFrame(\n",
    "    geometry=[ops.unary_union([p for p in world.loc['Netherlands', 'geometry'] if p.bounds[0] > 0])]\n",
    ").geometry.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "world.loc[['France'], 'geometry'] = gpd.GeoDataFrame(\n",
    "    geometry=[ops.unary_union([p for p in world.loc['France', 'geometry'] if p.bounds[1] > 0  and p.bounds[0] > -5])]\n",
    ").geometry.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "world = world.to_crs(epsg=3857)  # web mercator, so shapes are familiar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "countries = [\n",
    "    'Belgium', 'United Kingdom', 'New Zealand', 'Greece', 'Netherlands', 'Germany', 'Mexico', 'France', 'Peru',\n",
    "    'Singapore', 'Italy', 'Madagascar', 'Sri Lanka', 'Kenya', 'Turkey'\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And turn off the plots to make it easier to go out to file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using matplotlib backend: MacOSX\n"
     ]
    }
   ],
   "source": [
    "%matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "for country in sorted(world.index):\n",
    "\n",
    "    poly = world.loc[[country], 'geometry'].unary_union\n",
    "    \n",
    "    target_proportion = burned_area / country_areas.loc[country, 'area'] % 1\n",
    "    \n",
    "    _, miny, _, maxy = poly.bounds\n",
    "\n",
    "    solution = optimize.minimize_scalar(\n",
    "        minimise_area_proportional_difference,\n",
    "        method='Bounded',\n",
    "        bounds=(0, maxy - miny),\n",
    "        args=(target_proportion, poly)\n",
    "    )\n",
    "\n",
    "    solution.x\n",
    "\n",
    "    clipped_geom = clip_at_horizontal_line(poly, solution.x)\n",
    "\n",
    "    num_drawings = int(burned_area // country_areas.loc[country, 'area'] + 1)\n",
    "    num_cols = min(num_drawings, 10)  # at most 10 drawings wide\n",
    "    num_rows = int(math.ceil(num_drawings / 10))\n",
    "\n",
    "    fig, axes = plt.subplots(ncols=num_cols, nrows=num_rows, figsize=(10, 5))\n",
    "\n",
    "    images = [axes] if isinstance(axes, plt.Axes) else axes.ravel()\n",
    "    \n",
    "    for i, ax in zip(range(num_drawings - 1), images):\n",
    "        patch = PolygonPatch(\n",
    "            poly,\n",
    "            facecolor='red',\n",
    "            edgecolor='darkred',\n",
    "            linewidth=2\n",
    "        )\n",
    "        ax.add_patch(patch)\n",
    "        set_limits(ax, *poly.bounds)\n",
    "\n",
    "    ax = images[num_drawings - 1]\n",
    "\n",
    "\n",
    "    patch = PolygonPatch(\n",
    "        poly,\n",
    "        facecolor='lightgrey',\n",
    "        edgecolor='grey',\n",
    "        zorder=1\n",
    "    )\n",
    "    ax.add_patch(patch)\n",
    "\n",
    "    redpatch = PolygonPatch(\n",
    "        clipped_geom,\n",
    "        facecolor='red',\n",
    "        edgecolor='darkred',\n",
    "        linewidth=2,\n",
    "        zorder=2\n",
    "    )\n",
    "    ax.add_patch(redpatch)\n",
    "\n",
    "    set_limits(ax, *poly.bounds)\n",
    "    \n",
    "    for ax in images:\n",
    "        ax.axis('off')\n",
    "\n",
    "    fig.suptitle(\"Number of {}s burnt in Australia's 2019/2020 fires: {:.2f}\".format(\n",
    "        country, \n",
    "        burned_area / country_areas.loc[country, 'area']\n",
    "    ))\n",
    "    plt.savefig(f'country_fire_maps/{country}.png', dpi=600)\n",
    "    \n",
    "    plt.close()"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "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.7.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}